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ABSTRACT 

Fully developed turbulent and laminar flows through symmetric 
planar and axisymmetric expansions with heat transfer were modeled using 
a finite-difference discretization of the boundary- layer equations. By 
using the boundary- layer equations to model separated flow in place of 
the Navier-Stokes equations, computational effort was reduced permitting 
turbulence modeling studies to be economically carried out. The 
continuity and momentum equations were solved in a coupled manner. The 
validity of the once-through calculation scheme utilizing the FLARE 
approximation was studied by using a multiple sweep procedure in which 
the FLARE approximation is removed after the first sweep. 

For laminar constant property flow, the equations were 
nodimensionalized so that the solution was independent of Reynolds 
number. Two different dependent hydrodynamic variable sets were tried: 
the primitive variable set (u-v) , and the streamwise velocity stream 
function variable set (u-^>). The predictions of the boundary- layer 
equations were identical regardless of the variable set used. The 
predictions of the boundary- layer equations for parameters associated 
with the trapped eddy compared well with the predictions of the Navier- 
Stokes equations and experimental measurements for laminar isothermal 
flow 'd-en the Reynolds number was above 200 and the ratio of inlet to 
outlet channel diameter (width) was less than 1/3. The reattachment 
length and the flow field outside of the trapped eddy were well 
predicted for Reynolds numbers as low as twenty for laminar flow. 



xii 


The Boussinesq assumption was used to express the Reynolds stresses 
in terms of a turbulent viscosity. Near-wall algebraic turbulence 
models based on Prandtl 's-mixing- length model and the maximum Reynolds 
shear stress were compared. The near-wall models were used with the 
standard high -Reynolds -number k-e turbulence model. A low-turbulent- 
Reyno Ids -number k-e model was also investigated but found to be 
unsuitable for separated flow. The maximum-shear-stress near-wall model 
gave better predictions than the Prandtl -mixing- length models, 
especially for heat transfer. The predicted turbulent heat transfer is 
primarily dependent on the turbulence model used in the near-wall 
region. Globally iterating over the flow field had a more pronounced 
effect on the heat transfer solution than on the hydrodynamic solution. 
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I . INTRODUCTION 
A. Overview of Separated Flow 

Flow with separation is often encountered in fluid mechanics. 
Separation can occur due to adverse pressure gradients, as on the upper 
surface of an airfoil that has "stalled" or in a pipe with a sudden 
expansion or contraction. This latter type of separation occurs in 
engineering practice in heat exchangers and combustion chambers. The 
flow separation causes a region of flow reversal immediately downstream 
of the expansion that is sometimes referred to as a trapped eddy or 
trapped vortex due to its "swirling" nature. The region of flow 
reversal caused by the the flow separation can have an important effect 
on the flow field. 

Figure 1 introduces the physical nature of reattaching flow. The 
boundary layer at the step develops into a shear layer which has high 
levels of shearing stress. The recirculation region develops behind the 
face of the expansion. The streamline that divides the recirculating 
region from the rest of the flow is called the dividing streamline. The 
point where the dividing streamline meets the wall is called the point 
of reattachment. For "steady" turbulent flow, the point of reattachmsnt 
varies with time due to large turbulent eddies in the shear layer [1]. 
Thus, for turbulent flow, it may be more appropriate to define a 
reattachment region that extends 20 % of the reattachment length on both 
sides of the average distance to reattachment [2] . After reattachment. 
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sides of the average distance to reattachment [2]. After reattachment, 
a new sub-boundary layer begins to form. In the corner formed by the 
step and the wall downstream of the step, a smaller eddy is formed that 
rotates in a direction opposite to that of the larger eddy. An inviscid 
core above the shear layer may or may not exist depending on whether the 
flow is fully developed at the step. 

The present study deals with finite-difference solutions to partial 
differential equations that govern the velocity and temperature of fluid 
flow through rapid expansions. Since this geometry is one of the 
simpler geometries involving regions of flow reversal, it is of value as 
a test case to develop algorithms to be used with complex geometries. 

The Navier-Stokes equations are the general set of governing 
equations that are applicable for predicting the velocity field for this 
geometry. However, solving the Navier-Stokes equations requires more 
programming effort and computer time than solving equation sets that are 
more approximate. One such equation set consists of the boundary- layer 
equations . 

In general, the boundary- layer equations provide a good model of 
the flow field when the Reynolds number is very high. The boundary- 
layer equations can be obtained from the Navier-Stokes equations by 
assuming that the change of any variable of interest in the streamwise 
direction and the transverse velocity are both very small. The 
conservation of momentum equation in the transverse direction, y, 
reduces to a statement that the pressure varies only in the streamwise 
direction and is constant in the transverse direction (3p/3y=0). 
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For flows in which the boundary- layer assumptions are valid, the 
boundary- layer equations can be solved in place of the Navier-Stokes 
equations with relative ease. However, with rapid expansions, there 
exists a region of low Reynolds number flow immediately after the 
expansion in the recirculating eddy where the fluid velocity is low. 

The limitations of using the boundary-layer equations for this geometry 
has not been thoroughly addressed in previous studies. 

The boundary- layer equations are a parabolic set which means that 
no information from the downstream direction can influence the solution. 
Parabolic equation sets are applicable to flows in which a predominant 
direction of flow can be defined. This is in contrast to an elliptic 
equation set (such as the Navier-Stokes equations) in which events 
anywhere in the flow domain of interest can influence the solution. 
Elliptic equations can be used to model flow with no predominant flow 
direction. In regions of flow reversal, a predominant flow direction 
cannot be clearly defined. This indicates that flow regions with flow 
reversal are elliptic in nature. 

If the flow is not fully developed at the sudden expansion, a core 
region of inviscid fluid can be identified. Laplace’s equation 
(elliptic in type) can be effectively used to model flow in this 
inviscid core if the flow is irrotational , making it unnecessary to 
solve the full Navier-Stokes equations in this region. By using the 
boundary- layer equations in the viscous region near the wall and 
Laplace's equation in the inviscid free stream, elliptic effects can be 


included through the solution of Laplace's equation. The idea of using 
an equation set valid for the inviscid region and the boundary- layer 
equations near the wall is commonly known as viscous-inviscid 
interaction. However, if the flow is fully developed at the expansion, 
no internal core of inviscid fluid can be identified, so viscous- 
inviscid interaction cannot be used. For the flows predicted in this 
study, no inviscid core existed. 

When predicting turbulent flow, the governing equations are time 
averaged over a short period of time. As a result of averaging, extra 
terms are introduced into the momentum and energy equations and the 
original variables are replaced with time averaged variables. The extra 
terms in the momentum equations are called Reynolds stresses. 

Turbulence modeling of the Reynolds stresses in regions of flow reversal 
has been particularly challenging since many of the assumptions usually 
made in turbulence modeling are invalid when recirculation is present. 

The algebraic mixing length models neglect the diffusion and 
convection of turbulent kinetic energy and turbulent length scales. 

These algebraic models show good correlation wi"h experimental data for 
turbulence that is in equilibrium. In equilibrium boundary layers, the 
production of turbulent kinetic energy equals the dissipation of 
turbulent kinetic energy The turbulence encountered in recirculation 
flow is not equilibrium turbulence, so convection and diffusion of 
turbulence parameters need to be accounted for by using turbulence 
models more complex than the algebraic mixing length models. 
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More complex turbulence models that require the solution of 
partial-differential equations (PDEs) for the turbulence parameters are 
at the present time being refined. These models employ additional PDEs 
to take into account diffusion and convection of turbulence parameters. 
The turbulent kinetic energy equation (k-equation) is a PDE generally 
used 1' * „£ions of high turbulence Reynolds number. Some modifications 
ha’M been proposed to make it applicable for regions of low turbulence 
Refolds number and reversed flow. The k-equation requires the 
specification of a turbulent length scale either algebraically or 
through the solution of an ordinary differential equation (ODE) or 
another PDE. A popular PDE to provide this length scale is the 
e - equation (equation for the dissipation rate of turbulent kinetic 
energy). The e -equation, usually used in conjunction with the k- 
equation, is the weakest link of the k-e model [3]. The k-r. model is 
the most widely used two equation model of turbulence. The k-c model is 
used with the assumption that the Reynolds stresses are proportional to 
the mean rate of strain [4J . Other models do not make this assumption 
but impose PDEs for the Reynolds stresses themselves. These Reynolds 
stress models have not been widely used due to their complexity. 

However, algebraic equations that approximate the Reynolds stress PDE 
are often used in conjunction with the k-E model [ '] . 

If the temperature field is desired, the energy equation is used to 
solve for the temperature and to predict the heat transfer. If the flow 
is turbulent, the turbulent transport of energy must be modeled. As 
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with hydrodynamic turbulence models, there is a wide range of complexity 
involved with modeling the turbulent heat transport. In general, 
turbulence greatly enhances the heat transfer (and the skin friction). 

B. Literature Review 

In the last fifteen years, there has been a wealtn of new 
literature on the subject of reversed flow caused by a sudden 
contraction or expansion. However, for laminar flow heat transfer data 
sets are very scarce. There have been several studies of the heat 
transfer for turbulent flow, but until recently these studies have not 
included hydrodynamic data. Therefore, if differences appeared when 
comparing the temperature field or heat transfer rates with other data 
sets, one did not know if they were due to a faulty hydrodynamic 
solution or if they were due to problems with the solution of the energy 
equation itself. 

K Laminar rapid expansion studies 

a. Experimental hydrodynamic The laminar experimental 
hydrodynamic data sets are divided into two groups. Those dealing with 
axisymmetric and symmetric planar expansions are discussed in the 
symmetric section; those dealing with asymmetric planar expansions (also 
called backsteps or rearward-facing steps) are grouped in the asymmetric 
section 

1) Symmetric The first experimental study of an abrupt 


pipe expansion was carried out by Mactgno and Hung [5] who used flow 
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visualization with aluminum powder in oil. They considered only an 
expansion ratio of 1:2 (inlet diameter to outlet diameter) and found 
that the flow remained axisymmetric and stable for Reynolds numbers less 
than 4500. They also found that the distance from the step to the point 
of reattachment varied linearly with Reynolds number. 

Back and Roshko [6] used flow visualization to find the point of 
reattachment for flows in a circular channel abrupt expansion for 
Reynolds numbers between 20 and 4200. Dye was injected into water. 

They used a nearly uniform inlet profile with a very thin boundary layer 
at rhe expansion and reported that as the Reynolds number was increased, 
the distance to reattachment increased, reached a maximum, decreased, 
and then stayed constant when the flow became turbulent. The shear 
layer undulated for moderate Reynolds numbers with vortices appearing in 
it as the Reynolds number was increased further. For most of the range 
of Reynolds numbers studied, the flow was turbulent at reattachment. 

A 1:2 pipe expansion was studied by Iribarne et al. [7]. They 
found wave- like disturbances of the flow over a range of Reynolds 
numbers varying from 9 to 1355. The velocity profile at the step was 
nearly uniform with a very thin boundary layer. The shear layer 
downstream of the expansion developed wave-like disturbances. They 
reported that the frequency of the disturbances was a function of the 
square root of the Reynolds number. The irregular sinusoidal motion of 
the shear layer turned varicose for a Peynolds number of 350 which 
corresponded to the maximum reattachment length. There was considerable 
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turbulence at rcattachment when the Reynolds number was above 350. The 
mean flow field remained symmetric. 

Durst et al . [8], using a laser-Doppler anemometer, studied the 
flow through a symmetric 1:3 plane expansion over a range of Reynold 
numbers. The flow was fully developed at the inlet. Only for the 
lowest Reynolds number tried (56 based on the maximum velocity upstream 
of the expansion) was the flow symmetric. For Re of 252, the flow was 
stable but very asymmetric with additional separation regions downstream 
of the step. The two trapped eddies were of different length. An 
additional region of separation appeared after the shorter of the 
primary eddies for some Reynolds numbers. For Re above 252, the flow 
was unsteady. 

Feuerstein et al. [9] studied the fully developed flow through 
three different axisymmetric channel expansions using high speed 
photography of a fluid with tracer particles. They used the measured 
velocity profiles after reattachment as inlet conditions to a linearized 
boundary- layer problem to study the flow redevelopment. They gave a 
correlation for the distance to reattachment for d/D = 0.63 and 300 < Re 
< 1000 as 

£ /h = 0.048 Re*' 1 (1.1) 

r d 

Equation (1.1) predicts longer reattachment lengths than were predicted 


in this study. 
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Cherdron et al. [10] used a laser-Doppler anemometei to study flow 
that was fully developed at a 1:2 planar symmetric expansion. They 
found unequal reattachment lengths due to oscillations of the velocity 
if the Reynolds number based on the maximum inlet velocity was over 150. 
This caused antisymmetry due to the confinement of the flow. The 
reattachment length was close to an integral number of the wavelength of 
the disturbance in the shear layer. Disturbances at the lip were 
amplified in the shear layer to cause the unsteadiness. 

Restivo and Whitelaw [11] used the same channel as Durst et al. [8J 
but with a uniform inlet profile with a very thin boundary layer at the 
step. A laser-Doppler anemometer was used to measure the flow 
velocities. The flow was asymmetric for all Reynolds numbers studied 
(494 to 3865 based on the maximum inlet velocity). They studied a 1:2 
and a 1:3 symmetric planar expansion with a uniform velocity profile at 
the inlet. They reported a higher frequency of disturbance and so a 
shorter wavelength than reported by Cherdron et al. [10]. They believed 
that this caused the shorter distances to reattachment than those 
measured by Cherdron et al. [10] since the distance to reattachment 
remained the same integral number of wavelengths that Cherdron reported. 
However, Pollard [12] used the Navier-Stokes equations to study the 
effect of inlet conditions on flow over a step. He found that the 
reattachment length for a uniform inlet was less than that of a fully 
developed inlet for the same Reynolds number. Since Pollard's analysis 
did not take into account any wave- like disturbances, it is likely that 
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the shorter reattachment lengths measured by Restivo and Whitelaw [11] 
are only a result of the uniform inlet velocity profile and not the 
higher frequency of the shear layer disturbance. 

It is important to note that all of the symmetric laminar expansion 
experiments indicated that the flow has wave-like instabilities when the 
Reynolds number is of order 100 to 1000. The instabilities are worse 
for flows in which the boundary- layer thickness is approximately an 
order of magnitude less than the step height. Those having a boundary- 
layer thickness approximately the same order of magnitude as the step 
height are more stable [13]. The planar expansions are more unstable 
than the axisymmetric expansions because of the influence of the end 
walls of the planar expansions which introduce three-dimensional 
effects . 

When the boundary layer is fully developed at the expansion, 
experiments show that a laminar, two-dimensional, steady flow solution 
commonly exists. For a planar expansion the steady, two-dimensional 
flow was verified for a Reynolds number of 39 for a 1:3 expansion [8] 
and 100 for a 1:2 expansion [10]. For an axisymmetric expansion, 
steady, symmetric flow was observed for Reynolds numbers up to 4500 by 
Macagnc and Hung [ 5 ] . 

2) Asymmetric Since most of the asymmetric channel 
expansion (backward facing step) studies appearing in the literature 
involved a thin boundary layer at the step and an inviscid core or 
inviscid free stream, they were not used for comparison in the present 
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work. For this reason, some of the backward facing step data sets will 
only be briefly mentioned. 

Moore (14] was the first to study flow over a rearward -facing step. 
His study was carried out with air in a low speed wind tunnel. Leal and 
Acrivos [15] studied the effect of base bleed on the separated flow. 
O'Leary and Mueller [16] studied the flow over a backward facing step 
that had a thin boundary layer at the expansion using flow visualization 
with a water channel. Goldstein et al. [17] used a hot wire anemometer 
and raised some doubts concerning the earlier measurements of Moore 
[14]. The flow studied by Goldstein et al. had a thin boundary layer at 
the step. Denham and Patrick [18] used a laser anemometer to measure 
nominally fully developed flow over a back step. Matsui et al. [19] 
used flow visualization. Armaly and Durst [20], using a laser-Doppler 
anemometer to study the flow over a backward facing step, found 
additional regions of recirculation not found by Denham and Patrick 
[18]. Sinha et al. [21] studied flow in which the flow at the step had 
a thick boundary layer in relation to the step height but thin in 
relation to the channel height. 

Armaly et al. [22] found qualitative similarities between the 
stability of two-dimensional flow for symmetric and asymmetric 
expansions. They used a laser-Doppler anemometer with air as the fluid 
and found that the flow was two-dimensional for Reynolds numbers less 
than 400 and greater than 6600, corresponding to laminar and turbulent 
flow respectively. In the laminar range, the flow four step heights 
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upstream of the inlet was nearly parabolic. At the step, there was only 
a small variation from the parabolic profile. As the velocity was 
increased, the reattachment length increased when the flow was laminar, 
continued to increase as transition was reached, peaked, and then 
decreased until the flow was totally turbulent. This is similar to the 
findings of Back and Roshko [6] for an axisymmetric expansion. Once the 
flow was totally turbulent, the *-e attachment length remained constant. 

b. Experimental heat transfer It is important to note that 
there is no literature concerning experimental studies of heat transfer 
for symmetric laminar sudden expansions. Those investigators 
considering heat transfer for asymmetric planar expansions were Aung 
[23, 24] and Armaly et al. [25]. 

Aung [23,24] used a Mach-Zender interferometer to measure the heat 

transfer after a rearward-facing step. He used a uniform wall 

* 

temper ure and a fairly thick (6 /h ~ 1) boundary layer at separation. 
Aung gives reattachment lengths, free stream velocities, heat transfer 
coefficients, average Stanton numbers, temperature distributions, and 
the following correlation for the average Stanton number. 

St = 0.787 Re"' 55 (s/x ) 0,72 (1.2) 

s s 

where Re g is the Reynolds number based on the step height, s, and the 
free stream velocity above the step, and x g is the distance from the 
step. Aung found chat the local heat transfer coefficient increased 
: 'foi.-riically in the separated region but was always lower than that for 
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flow over a flat plate. Regardless of this, for transitional flow, the 
total heat transfer can be significantly higher than for flow over a 
flat plate. The maximum heat transfer was downstream of reattachment. 
Aung also found that the distance to reattachment peaked while the flow 
was in transition, then decreased until it became turbulent. The stream 
line curvature before the step was said to have increased the heat 
transfer upstream of the step. 

Armaly et al. (25] measured the momentum, heat and mass transfer in 
backward- facing step flows that were fully developed at the step. Their 
reported heat transfer data show that the Reynolds analogy is invalid in 
the region of reversed flow. At Reynolds numbers near 400, the laminar 
flow became three-dimensional. 

c. Numerical hydrodynamic The numerous numerical solutions for 
flows through rapid expansions can best be summarized in tabular form. 
Table 1 shows the majority of the available computational solutions for 
laminar flows through rapid expansions to date. 

1) Symmetric Hung [26] was the first to calculate laminar 
two-dimensional flows of this type using the Navier-Stokes equations. 

He predicted the flow through a 1:2 symmetric planar expansion and a 1:2 
axisymmetric expansion for Reynolds numbers less than 360. He found 
that the distance from the step to the point of reattachment varied 
linearly with Reynolds number. Macagno and Hung [5] compared the 
axisymmetric Navier-Stokes solutions with their own measurements. The 
correlation was excellent. Hung [26] and Macagno and Hung [5] used the 
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Table 1. Summary of laminar numerical solutions for sudden expansions 


Author (s) 

Method 

Expansion Type 

Inlet Condition 

Acrivos & 
Schrader [27] 

BL, Finite Dif. 

Sym. Planar 

Parabolic, 

Uniform 

Agarwal [28] 

NS, Finite Dif. 

Sym. Planar 

Parabolic 

Armaly et al. 
[22] 

NS, Finite Dif. 

Back-Step 

Parabolic 

Atkins et al. 
[29] 

NS, Finite Dif. 

Back-Step 

Experimental 3 

Parabolic 

Chen et al. 
[30] 

NS, Finite 
Analytic 

Back-Step 

Parabolic 

Durst et al. 
[8] 

NS, Finite Dif. 

Sym. Planar 

Parabolic 

Giaquinta [31] 

NS, Finite Dif. 

Sym. Planar 

Uniform 

Hackman et al . 
[32] 

NS, Finite Volume 

Back-Step 

Parabolic, 

Experimental 

Halim & Hafez 
[33] 

PPNS, Finite Dif. 

Sym. Planar 

Parabolic 

Hall & Pletcher 

Viscous/ Inviscid 

Back-Step 

Thin boundary 

[34] 

Interaction 


layer 

Hung, T.K. [26] 

NS, Finite Dif. 

Sym. Planar, 
Axisymmetric 

Parabolic 

Hutton & Smith 
[35] 

NS, Finite Elem. 

Back-Step 

Parabolic 

Huyakorn et al. 
[36] 

NS, Finite Elem. 

Back-Step 

Parabolic 

Kumar & Yajnik 
[37] 

BL, eigenfunc- 
tion expansion 

Sym. Planar 

Parabolic 

Kwon et al. 
[38] 

BL, Finite Dif. 

Sym. Planar 

Parabolic, 

Experimental 

Kwon & Pletcher 

Viscous/ Inviscid 

Sym. Planar 

Par at j lie, Thin 

[39] 

Interaction 


boundary layer, 
Experimental 

Leschziner 

[40] 

NS, Finite Dif. 

Back-Step , 
Axisymmetric 

Parabolic 

Macagno & Hung 
r s ] 

NS, Finite Dif. 

Axisymmetric 

Parabolic 

l 1 J 

Morihara [41] 

NS, Finite Dif. 

Sym. Planar 

Parabolic 

Oosthuizen 

[42] 

BL, Finite Dif. 

Back-Step 

Parabolic, 

4th order poly 

Osswald et al. 
[43] 

NS, Finite Dif. 

Sym. Planar 
Axisymmetric 

Parabolic 


All the inlet conditions listed as experimental where nearly 
parabolic . 
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Table 1. (continued) 


Author (s) 

Method 

Expansion Type 

Inlet Condition 

Plotkin [44] 

BL, Fourier 
Series approx. 

Sym. Planar 

Parabolic 

Pollard [12,45] 

NS, Finite Dif. 

Axisymmetric 

Parabolic 

Uniform 

Roache & 
Mueller [46] 

NS, Finite Dif. 

Back-Step 

4th order 
Polhausen poly 

Taylor et al. 
[47] 

NS, Finite Elem. 

Back-Step 

Parabolic 

Experimental 

Thomas et al . 
[48] 

NS, Finite Elem. 

Back-Step 

Parabolic 


stream function and the vorticity as variables in both a steady 

and an unsteady approach. 

Morihara [41] was next to solve the Navier-Stokes equations for 
flow through symmetric channels. Instead of using the vorticity-stream 
function variables he developed a primitive variable technique with the 
pressure terms eliminated from the equations. He coupled the continuity 
eqi’.ation with the momentum equations and thus circumvented the need to 
use a relaxation parameter. The algorithm appeared stable but required 
large computer storage if the coefficient matrix was not broken into 
smaller subraatrices. He predicted shorter reattachment lengths than 
those predicted by Hung [26] for a 1:2 symmetric planar expansion. 

Durst et al. [8] solved the Navier-Stokes equation^ using the 
vorticity-stream function ..ethod first developed by Gosman et. al . [49] 
to compare with their experimental measurements through a 1:1 (d:D as in 
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Fig. 2) symmetric planar expansion. However, since the flow was 
symmetric and two-dimensional for only the lowest Reynolds number tried 
(56 based on the maximum inlet velocity), this was the only comparison 
they were able to make. 

Giaquinta (31) used two types of differencing molecules to solve 
the unsteady ip-u form of the Navier-Stokes equations as a model for the 
flow through a 2:5 symmetric planar expansion. He studied the 
difference between an explicit time method and an implicit time method. 
The flow was started from rest and allowed to reach a steady state. He 
found that the implicit method was good for long time analyses, 
particularly after flow initiation. The explicit time method was good 
for predicting sudden changes in the fluid motion such as at start up. 
The inlet velocity profile was uniform. The flows studied had Reynolds 
numbers of 10 and 100. 

Leschziner [40] used primitive variables to test the predictions of 
the Navier-Stokes equations for three different types of finite 
differencing of the convective terms for an axisymmetric and asymmetric 
planar expansion. His predicted reattachment length compared well with 
the measurements of Macagno and Hung [5] for an axisymmetric expansion. 
He stated that artificial diffusion due to skewness of grid and 
streamlines is unimportant for laminar flow. 

Pollard [12,45] used a primitive variable formulation of the 
Navier-Stokes equations to predict the flow through axisymmetric 
expansions by finite-difference discretization. His computational 
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algorithm was similar to the SIMPLE algorithm developed by Patankar and 
Spalding [50]. He studied the effect of varying the Reynolds number, 
expansion ratio, and inlet profile and found that the distance to 
reattachment varied linearly with Reynolds number and nonlinearly with 
expansion ratio. Reattachment lengths were shorter for uniform inlet 
profiles than for parabolic inlet profiles. It should be noted that for 
some Reynolds numbers his predictions of c^ exceeded the known fully 
developed c^ values by as much as 11%. His predicted reattachment 
lengths compared well with those predicted by Macagno and Hung [5], 

Agarwal [28] used a third-order accurate upwind differencing scheme 
to solve the tfi-w form of the Nsvier-Stokes equations. However, his 
predictions do not compare well with the measurements of Durst et al. 

[8] for a 1:3 symmetric planar expansion. For a 1:2 symmetric planar 
expansion, the results compared well with the predictions of Kumar and 
Yajnik [37] and Kwon and Pletcher [39]. The algorithm was stable and 
accurate for high Reynolds numbers. 

Osswald et al. [43] used a direct implicit time-dependent technique 
to solve the vorticity-stream function form of the Navier-Stokes 
equations to predict the flow through a 1:3 symmetric planar expansion 
and a 1:2 axisymmetric expansion. A generalized orthogonal coordinate 
system with a cluster conformal transformation technique packed the grid 
in the regions where the length scale was shorter. An ADI (alternating 
direction implicit) scheme was used to solve the vorticity transport 
equation; a block-Gaussian elimination scheme was used to solve the 


stream function equation. 
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Those who used the partially-parabolized Navier-Stokes equations 
(PPNS) include Madavan [51] and Chiu [52] . They neglected the 
streamwise diffusion terms in the Navier-Stokes equations. 

Only recently have the boundary- layer equations been used to 
predict the flow through the large regions of reversed flow behind 
sudden expansions. Acrivos and Schrader [27] felt that the boundary- 
layer equations were not valid near the step. In this "near region" 
close to the s_ep, they stated that the flow was inviscid in nature, the 
length of the inviscid region being of the same order of magnitude as 
the inverse of the Reynolds number. To overcome this, they modified the 
initial velocity profile at the step to account for upstream influence 
in the near region. Acrivos and Schrader used the unsteady boundary- 
layer equations by marching in time until a steady state solution was 
reached. They added damping to the viscous term in the momentum 
equation to suppress instabilities in the. finite-difference procedure. 
They predicted the flow through symmetric planar expansions for various 
expansion ratios. 

Kumar and Yajnik [37] reported that by properly scaling the 
coordinate parallel to the channel centerline, the region of inviscid 
flow was shrunk to zero. They argued that a form of the boundary- layer 
equations, as a set of limit equations to the Navier-Stokes equations, 
is applicable for this geometry for sufficiently high Reynolds numbers. 
Kumar and Yajnik used an eigenfunction expansion to reduce the set of 
partial differential equations to a set of ordinary differential 
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equations. They predicted the flow through symmetric planar expansions. 
Plotkin [44] solved the boundary- layer equations in a way similar to 
Kumar and Yajnik [37]. However, instead of an eigenfunction expansion, 
Plotkin used a concise Fourier series approximation technique. 

Plotkin' s results were similar to those of Kumar and Yajnik. Both 
Plotkin [44] and Kumar and Yajnik [37] reported singular behavior when 
attempting to predict separated flow if too many expansion terms were 
included. 

Kwon et al. [38] solved the boundary- layer equations with a once 
through marching procedure. In regions of reversed flow, the FLARE 
approximation [ 53 ] removed the streamwise convective terms. The 
momentum and continuity equations were solved in coupled manner. The 
predicted reattachment length compared well with that predicted by the 
Navier-Stokes equations and that measured experimentally. 

2) Asymmetric Roache and Mueller [46] used a finite- 
difference discretization of the unsteady Navier-Stokes equations to 
predict the flow field passing over an asymmetric planar expansion or 
back-step. They used a ip-u form of the Navier-Stokes equations. They 
marched explicitly in time until the solution stabilized. 

Atkins et al. [29] solved the ip-u form of the Navier-Stokes 
equations with finite differences. Their laminar predictions used the 
experimer tal, inlet profile of Denham [54] (2:3 asymmetric expansion). 
They tried upwind and central differencing of the convective terms. The 
predictions of both differencing schemes were very close to the 
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aeasurements of Denham. For Reynolds numbers lower than those measured 
by Denham, the upwind differences predicted reattachment lengths and 
eddy intensities that were 8% less than those predicted by the central 
differences . 

Leschziner [40] used primitive variables to test three methods of 
finite differencing the convective terms of the Navier-Stokes equations. 
He compared his predictions with the measurements of Denham and Patrick 
[18] for 2:3 asymmetric expansion. He predicted the minimum stream 
function value in the eddy measured by Denham and Patrick but over 
predicted the measured reattachment length. He stated that artificial 
diffusion due to skewness of grid and streamlines is unimportant for 
laminar flow. 

The TEACH code [55] was used by Armaly et al. [22] to solve the 
Navier-Stokes equations. Their measurements and predictions compared 
well as long as they used their measured inlet condition in their 
computations and as long as the experimental flow remained two- 
dimensional. The back-step they studied had an exparsion ratio of 
1:1.94. 

Hackman et al. [32] used a finite volume discretization of the 
primitive variable Navier-Stokes equations to test two different 
differencing schemes. They predicted the flow over a backward facing 
step and compared the results with Denham and Patrick [18] for Reynolds 
numbers of 73 and 229 based on tne step height and average inlet 
velocity. When the measured inlet profile at the step was used as the 
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inlet condition for computations, the numerical predictions compared 
well with experimental measurements. Their predictions did not compare 
well with measurements when a parabolic inlet condition was used. 

There have been three studies that used viscous-inviscid 
interaction for back-step flows that had a thin boundary layer at the 
step. Kwon and Pletcher [39] developed a hydrodynamic solution method. 
In the viscous region, they used the boundary- la>er equations with the 
FLARE approximation ' J] ; in the inviscid region they solved Laplace's 
equation. Halim and Hafez [33] solved a fourth-order equation for the 
stream function in the viscous region derived from the PPNS equations. 
Halim and Hafez introduced an implicit coupling procedure for coupling 
the viscous and inviscid solutions. 

There have been several finite element solutions of the Navier- 
Stokes equations for flow over back-steps. The usual method of using 
the Galerkin formulation with weighted residuals that is so widely US ed 
in structural mechanics often produces oscillations in the solution. 

This is because the convection terms cannot be easily treated by the 
symmetrical operators that are commonly used in structural mechanics. 

Huyakorn et al. [36] predicted the fully developed flow through an 
asymmetric expansion using different interpolation elements. Hutton and 
Smith [35] predicted the same flow case. 

Taylor et al. [47] used a weighted residual finite element 
discretization of the primitive variable Navier-Stokes equations to 
predict the laminar and turbulent flow over a back-step. They compared 
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their predictions with the experimental measurements of Denham and 
Patrick [18] and the finite-difference predictions of Atkins et al. 

[29]. The predictions of Taylor et al . did not match the measurements 
as well as the finite-difference predictions of Atkins et al. 

Thomas et al. [48] predicted the laminar and turbulent flow over a 
back-step using primitive variables. They had to substitute upwind 
weighting functions for the Galerkin weighting functions to get 
convergence. This is similar to the need to sometimes use upwind 
differencing instead of central differencing of the Navier-Stokes 
equations to reach a stable finite-difference solution. The variable 
coefficients of the convective terms were set to the previous global 
iteration values. This provided a means of linearization and a means of 
global iteration. The predictions of Thomas et al. compare reasonably 
well with the measurements of Denham and Patrick [lfij and the finite- 
difference predictions of Atkins et al . [29]. 

Chen et al. [30] used what they termed a finite-analytic procedure 
to solve the Navier-Stokes equations. They predicted the flow over a 
2:3 back-step. They predicted more mass trapped in the eddy behind the 
step than that measured by Denham and Patrick [18]. 

Chiu [52] used the PPNS equations to predict the incompressible 
flow over a back-step (see Numerical-symmetric section). His 
predictions are in good agreement with the predictions of the full 
Navier-Stokes equations and experimental measurements, indicating that 
the strearawise diffusion terms are unimportant for moderate Reynolds 


numbers . 
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Kwon and Pletcher [39] and Oosthuizen [42] used a finite-difference 
formulation of the boundary- layer equations to predict fully developed 
flow over an asymmetric expansion. Oosthuizen was able to predict t'»e 
additional region of separation on the wall opposite the step that was 
first reported by Armaly and Durst [20]. Both Kwon and Pletcher [39] 
and Oosthuizen [42] used the FLARE approximation to march through 
regions of reversed flow. 

The previous studies that used the boundary- layer equations have 
provided only isolated comparisons indicating that the solution to the 
boundary- layer equations may provide useful information for some sudden 
expansion flows. However, the limitations of the boundary- layer 
equations have not been clearly defined. 

d. Numerical heat transfer Hall and Pletcher [34] modified the 
algorithm of Kwon and Pletcher [39] to include a solution of the energy 
equation. Theirs is the only laminar solution of the energy equation 
for flow over a back-step in the literature. 

2. Turbulent rapid expansion studies 

The experimental studies are discussed first, followed by the 
numerical studies. Both the experiment* 1 and numerical studies are 
divided into two groups according to whether or not they include heat 
transfer data. 

Turbulent measurements are more difficult to make than laminar flow 
measurements. The early studies used flow visualization, pitot tubes, 
hot wires, and surface pressure transducers to measure flow parameters. 
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At first, these data provided mostly qualitative, but important, flow 
field information. Later as hot wire techniques were refined, more 
accurate measurements were taken. With the introduction of laser- 
Doppler anemometers and pulsed -wire anemometers, the measurements again 
increased in accuracy. At the present time, refinements are still being 
made to increase the accuracy of turbulence measurements [2j. 

a. Experimental hydrodynamic The experimental -hydrodynamic 
studies are divided into two groups. Those dealing with flow through an 
axisymmetric or a planar symmetric channel expansion are grouped in the 
symmetric category; those dealing with flow over a backward- facing step 
are grouped in the asymmetric category. 

1) Symmetric The studies giving only hydrodynamic data 
for symmetric sudden expansions are summarized as follows: 

Flow visualization: Drewry [56] 

Pitot tubes: Kangovi and Page [57], Ha Mini J . Chassaing [58], 

Mehta [59] 

Hot-wire anemometers: Abbot and Kline [60], Chaturvedi [61], Ha Minh and 
Chassaing [58], Mehta [59] 

Laser-Doppler anemometers: Moon and Rudinger [62], Freeman [63], 

Smyth [64], Lu [65], Stevenson et al. [66] 

Abbot and Kline [60], Smyth [64], and Mehta [59] all studied the 
flow through symmetric (aouble-sided) planar expansions. Abbot and 
Kline [60] were the first to experimentally study the flow in asymmetric 
and symmetric planar expansions using hot wire probes. Smyth [64] was 
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the first to study the flow through a symmetric planar expansion using a 
laser-Doppler anemometer. He provided mean velocity, turbulent kinetic 
energy, Reynolds stress, and fluctuating velocity profiles. The flow 
was fully developed and turbulent at the step. He found no appreciable 
flow asymmetries. Mehta [59] found asymmetric unsteady flow through a 
symmetric expansion when d/D was greater than 1.5. For d/D smaller than 
1.5, the flew remained steady and symmetric. The asymmetries may have 
been due to three dimensional effects caused by the small channel aspect 
ratio (the ratio of the channel height to the width was only 1/4). 

Mehta used hot wires and pitot tubes . 

There have been numerous axisymmetric rapid expansion studies. 
Chaturvedi [61] used pitot tubes and hot wire anemometry to study 
axisymmetric expansion flow. He provided velocity, pressure, and 
turbulence data for different step face wall angles. 

Moon and Rudinger [62] studied fully developed turbulent flow 

through a axisymmetric expansion using a laser-Doppier anemometer (LDA) . 

They published velocity profiles, centerline velocity curves, and eddy 
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shape diagrams. For the range of Reynolds numbers studied (10 -10 
based on the inlet diameter), the reattachment length was approximately 
1.25 outlet tube diameters from the step and independent of Reynolds 
number. They also compared their measurements with numerical 
predictions . 

Freeman [6.3] studied flows that had hot and cold co-flowing streams 
at the expansion. He measured the reattachment length as 5 step heights 
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from the expansion for Reynolds numbers between 20,000 and 40,000 using 
a laser anemometer. For this range of Reynolds numbers, he found that 
the temperature profiles were independent of the Reynolds number. 

Kangovi and Page [57] used pitot tubes to measure the flow through 
an axisymmetric sudden expansion. They measured a reattachment length 
of 8 step heights. 

Ha Minh and Chassaing [58] used hot wires and pitot tubes to 
measure flow that had a very thin boundary layer through a 1:2 pipe 
expansion. They reported a turbulence intensity that was 19% of the 
centerline velocity in the reattachment region. The turbulence 
intensity decayed rapidly downstream of reattachment. They measured 
reattachment at nine step heights from the step. 

Stevenson et al. [66] measured the velocities of flow through an 
axisymmetric expansion and an asymmetrical backstep using a laser 
anemometer. They used frequency shifting and control of the seeding 
density to eliminate bias errors. The axisymmetric case had a short 
entry length before the expansion. The axisymmetric inlet profile was 
measured with a pitot tube and found to be "very flat". They found the 
peak Reynolds stress was at the edge of the recirculation zone. They 
published turbulent kinetic energy, Reynolds stress, and some 
fluctuating velocity data. They also predicted the flow using a 
modified version of the SIMPLE computer code with a k-e turbulence 
model. The strength of their paper is in the experimental measurements. 
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Drewry [56J and Lu [65] studied flows having a constriction at the 
inlet of a relatively short pipe. Lu [65] used a laser anemometer to 
determine the velocities in a short tube with constrictions at both 
ends. He provided the centerline velocity and a few velocity profiles. 
The inlet conditions, which were not provided, were not that of fully 
developed turbulent flow. Drewry used flow visualization. 

2) Asymmetric Eaton and Johnston [1] gave an excellent 
summary of the backward facing step experiments that were done before 
1931. Rather than repeat their work, the main points of their summary 
will now be listed and a detailed sketch of the experiments they 
summarized will not be given. 

Upon reviewing the early experiments, Eaton and Johnston made the 
following conclusions. 

• The reattaching shear layer is like a free shear layer except 
on the :, wall side" where the flow is highly turbulent. Since 
the outer part of the reattaching shear layer is similar to a 
free plane shear layer, turbulence models that work well for 
free shear layers should work well for the reattaching shear 
layer except near the wall and reattachment point. Some 
refinement will undoubtedly have to be done for these regions. 

The outer part of the reattaching layer retains the 
characteristics of a free shear layer for as many as 50 step 
heights downstream of the step. 
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• There is a rapid decrease in Reynolds stresses near 
reattachment that is caused by either streamline curvature, 
pressure gradient, wall interaction or some combination of the 
above . 

• The maximum reversed flow velocity is usually over 20% of the 
free stream velocity. 

• The boundary layer state (turbulent or not) at the step has an 
important effect on the downstream flow. For fully turbulent 
boundary layers, the flow is independent of the Reynolds 
number. However, it is not clear how much the boundary layer 
thickness affects the flow. 

• Increasing the free stream turbulence tends to decrease 
leattachment lengths; increasing the expansion ratio tends to 
increase reattachment lengths. 

• Hot-wire probes tend to measure lower turbulence quantities 
than laser anemometers. It was concluded that the hot-wire 
anemometers under measure turbulence quantities. 

• The y-position of the maximum turbulence intensity moves toward 
the wall as reattachment is approached; it moves away from the 
wall as reactachmert is passed. 

• There still remains controversy on whether spanwise vortices 
are the dominant structure in the plane mixing layer. Eaton 
and Johnston concluded that this is the case. They concluded 
that these vortices are the cause of the unsteady reattachment 
point noted in several studies. 
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• There is also controversy on what happens to the large eddies 
at reattachment . Some feel these eddies are ripped apart while 
others feel that some of the large eddies go upstream and some 
go downstream. 

Since Eaton and Johnston's review of backward facing step 
experiments, there have been several new studies that have added insight 
to the physics of the flow. Since the review. Driver and Seegmiller 
[67], Pronchick and Kline [68], Adams [69], and Stevenson et al. [66] 
have used laser anemometers to study the recirculating flow. Pronchick 
and Kline also used flow visualization. Cheun et al. [70] and Mess and 
Baker [71] used pulsed wire anemometers; Westphal et al. [2] used a 
pulsed wall probe. 

Along with using a pulsed wire anemometer, Moss and Baker [71] also 
measured the surface pressure for flow over small protuberances and a 
backward facing step in a large wind tunnel. They found that "the line 
of peak stresses diverges progressively outwards from the dividing 
streamline with values rising to a maximum before falling away with 
reattachment", which is in agreement with the summary of Eaton and 
Johnston [1]. 

Cheun et al. [701 studied the effects of the free stream turbulence 
and the boundary layer thickness at the step. They reported that for 
their experiments, the free stream turbulence had little effect on the 
flow. This is different from the conclusion reached by Eaton and 
Johnston [1] in their review. Cheun et al. also reported that the 
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thinner the boundary layer at the step, the shorter the reattachment 
length. 

Chandrsuda and Bradshaw [72] studied the turbulent stresses and 
turbulent energy balance in the reattachment region with hot-wire 
anemometers. They concluded, in agreement with everyone else, that the 
shear layer upstream of reattachment is similar to a plane mixing layer. 
They said that the change of the turbulence structure near reattachment 
is due to the confinement of the large eddies. They felt that an 
accurate turbulence model should have a fairly sophisticated model for 
the triple products, i.e., a triple product transport equation and that 
the dissipation equation should include a wall effect term. 

Driver and Seegmiller [67] measured the thickness of an oil film 
with a laser to determine c^. They reported a sudden increase in 
turbulent stress after separation which started to decrease two step 
heights before reattachment. The triple cross products rapidly 
disappeared at reattachment which suggested that the eddies were being 
torn apart. They reported that the wall side of the shear layer was 
highly turbulent. In comparing with numerical predictions, they 
concluded that using an algebraic stress model improved predictions 
because streamline curvature was important. 

Pronchick and Kline [68] hopefully settled the dispute concerning 
what happens to the eddies at reattachment . According to Pronchick and 
Kline, the two earlier theories concerning the fate of the large eddies 
are both correct. Pronchick and Kline concluded that some eddies are 
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broken apart it reattachment and some are either swept upstream or 
downstream. Those swept downstream are the cause of the slow recovery 
of the flow to typical boundary layer flow. They reported that the 
unsteadiness of the reattachment point was caused by three-dimensional 
eddies . 

Westphal et al . [2] mainly used a newly developed pulsed wall probe 
to study the flow in the back-flow region. They found a strong 
dependence of the rcattachment length on the boundary layer thickness . 
The flous they studied became similar when the x-coordinate was scaled 
about the reattachment point. They felt that there was a laminar- like 
thin region of strong backflow next to the surface upstream of 
reattachment that was not similar to attached turbulent boundary layers. 

The investigation by Stevenson et al. (66] was mainly to study the 
errors in velocity biasing when using laser anemometers for turbulen' 
flow. Velocity bias is due to the fact that for turbulent flow, more 
particles per unit time traverse the probe volume when the velocity is 
high than when it is low. Velocity bias can be overcome by high uniform 
seeding density and equal time velocity sampling. They found the peak 
Reynolds stress at the edge of the recirculation zone. The channel was 
so narrow that the measurements were not truly two-dimensional since the 
side walls undoubtedly affected the flow. For this reason, the 
asymmetric planar expansion measurements cannot be used as a comparison 
for two-dimensional numerical predictions. 
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Walterick et al. [73] found relatively high free stream turbulence 
in their backward facing step measurements. They concluded this was 
caused by the unsteady motion of the shear layer that was noticed at the 
reatcachment point. However, it may have been due to inadequate flow 
conditioning upstream of the test section [74]. The importance of free 
stream turbulence is still unclear [1,70]. 

Lamb and McCotter [75] made surface pressure measurements over a 
small step and protuberances in a large wind tunnel. They were able to 
correlate the pressure in the recirculation region using reference 
values at reattachment and the point of minimum pressure. 

Adams [69] found that a laminar boundary layer at separation gave 
shorter reattachment lengths than a turbulent boundary layer. The flow 
remained Reynolds number independent for the range of Reynolds numbers 
studied (under 36,000 based on the step height). He, like Westphal et 
al . [2], argued that the boundary layer in the recirculating region was 
laminar-like. He found no bursting mechanism like that observed i^ 
typical turbulent boundary layers. 

The reversed flow studied by Simpson et al. [76] and Simpson [77] 
was for flow that had separation induced by a pressure gradient on a 
flat plate but has application to rapid expansion flows. Simpson et al. 
[76] reported that the eddy viscosity and mixing length models are poor 
in the separation region. They claim that the normal law of the wall in 
u + and y + coordinates is not valid in the recirculation bubble near the 
wall. They say that the turbulence stress must be modelled according to 
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the turbulence structure and not by the mean velocity gradients. They 
state that the reason the eddy viscosity models are not good in the 
reversed flow is because 3u/9y is based on averaging the large turbulent 
fluctuations. These averages are not m eaningful due to the relative 
magnitudes of the instantaneous velocity and the time averaged velocity. 
Simpson [77] divides the flow into three layers: (1) a viscous layer 
near the wall, (2) an intermediate layer that acts as an overlap or 
buffer, and (3) an outer layer which is part of the large scale outer 
flow. For the viscous layer near the wall, he proposes the following 
equation 

0-3) 

N 

where N is the location of the maximum negative velocity in the bubble 
and u^ is the absolute value of the maximum negative velocity. Using 
Eq. (1.3), he developed wall functions that can be used in place of the 
law of the wall. 

b. Exper i mental heat transfer There have been several reviews 
of the heat transfer data for turbulent separated flows. Hanson and 
Richardson [78] and Chilcott [79] were the earliest. Fletcher et al. 

[80] reviewed a large number of papers published previous to 1974 for 
both subsonic and supersonic separating reattaching flows for various 
geometries. Aung and Watkins [81] and Aung [24] reviewed the turbulent 
subsonic heat transfer studies in 1978 and 1983 respectively. 
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Aung's review [24] emphasized the studies of recirculating flows 
over steps and cavities for laminar, transitional, and turbulent flows 
that were done after 1978. He noted that downstream of the step, the 
pressure reaches a low value and then remains almost constant in the 
axial direction for some distance. He concluded that the curvature 
effect of the streamlines enhances the heat transfer. The studies he 
reviewed indicated that the thermal boundary layer redeveloped to a form 
that was typical of flat plate flow quicker than the velocity field 
redeveloped, ’^'e also pointed out the difference between laminar and 
turbulent separation heat transfer. For laminar heat transfer, the heat 
transfer rate in the recirculation bubble is always less that that for a 
flat plate. It starts low and increases monotonically . On the other 
hand, turbulent heat transfer greatly exceeds that of fully developed 
flat plate flow in the recirculating region. The turbulent heat 
transfer rate peaks at reattachment then drops to fully developed flat 
plate or channel flow values. Aung states that high levels of 
turbulence are generated in the shear layer where the turbulence 
dissipation is small. The dissipation remains small until the length 
scale decreases due to the flow approaching the wall. This is why the 
turbulent stresses increase in the reattaching shear layer only to 
decrease when reattachment is approached. 

The different heat transfer studies listed by geometry are: 

• Axisymraetrio expansion: Ede et al. [82], Krall and Sparrow 
[83], Zemanick and Dougalt [84], Runchal [85], Back et al. 

[86], Kang et al. [87], Sparrow and O'Brien [88], 
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Amano et al. [89], and Baughn et al . [90] 

• Planar symmetric: Filetti and Kays [91], Seki et al. [92], and 
Seki et al. [93] 

• Asymmetric rearward- facing step: Seban et al. [94], Seban [95], 
Aung and Goldstein [96], Kottke *97], and Vogel and Eaton [74] 

There is some disagreement concerning the effect of Reynolds number 

on the heat transfer. Most of the planar symmetric and exisymmetric 

expansion studies indicated that the maximum Nusselt number is 
2/3 

proportional to Re [83, 84, 90, 92, 93]. However, Amano et al. [89] 
found that for the smallest expansion ratio they studied (d/D = 0.195), 

the Nusselt number was not a function of the Reynolds number and varied 

0 5 ^ 

according to Re ' for the larger two expansion ratios that they 

studied. Ede et al. [82], who used water as a fluid with heatin'* 

upstream of the expansion as well as downstream, found the variation of 

the convective heat transfer coefficient substantially independent of Re 

over a wide range of Reynolds numbers (3700-45,000). Filetti and Kays 

[91] predicted two different reattachment lengths for their symmetric 

planar expansion and so measured the Nusselt number as being 

proportional to Re” 1 , wnere the measured values of m bracketed 2/3. 

Filetti and Kays [91] reported that the Nusselt number was proportional 

to Re°'*\ For rearward- facing step flows, Seban et al. [94] found the 

0 8 

maximum Nusselt number proportional to Re . Sparrow and O'Brien [88] 
measured the heat transfer along the face of an axisymmetric expansion 
and concluded that for high Reynolds numbers, the heat transfer rate 
became less and less dependent on Re. 
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Most of the studies indicate that the maximum Nusselt number is at 
the reattachment point [83, 93, 94]. Vogel and Eaton [74] used a newly 
developed pulsed wall probe that measures to verify that the point of 
maximum heat transfer and zero shear stress correspond. However, Kang 
et al. [87] using flow visualization measured it to be upstream of 
reattachraent by 15% of the reattachment length. The value of the 
maximum Nusselt number was higher than the fully developed Nusselt 
number by a factor of 2.7 to 11, depending on the reference. The 
maximum Nussult number is very geometry dependent [84] . 

Zemanick and Dougall [84] reported a small effect of Re on the 
reattachment point except for very high and very low Reynolds numbers. 
For the very high Reynolds numbers, compressibility might have had an 
effect. Tney found that reattachment was a function of the expansion 
ratio with smaller values of d/D giving longer reattachment lengths. 
Baughn et al. [90], whose test was very similar to that of Zemanick and 
Dougall, found that for a given expansion size, the bubble length varied 
little with Re. Both the studies by Zemanick and Dougall and Baughn et 
al. utilized fully developed flow at the axisymmetric expansion. Krall 
and Sparrow [83] measured a reattachment length of 1.25 to 2.5 diameters 
from the step. The reason why some studies showed little effect of 
Reynolds number and expansion size on the distance to maximum Nusselt 
number and others showed a large effect is probably due to differences 
in the boundary layer thickness at the expansion. 
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Several investigators found evidence of a small eddy in the corner 
near the step rotating in a direction opposite to that of the larger 
eddy. Zemanick and Dougall (84] and Baughn et al. [90] both conclude 
that since there is a local maximum in the heat transfer rate very near 
the step, this indicates the existence of the small eddy. Sparrow and 
O'Brien [88] said their napthalene sublimation method made the existence 
of a secondary eddy obvious . 

The studies point out that the heat transfer through the near wall 
layer does not behave like that through a normal turbulent boundary 
layer in the reattachment and redevelopment regions. Seban [95] found 
large temperature gradients near the wall downstream of the step. Aung 
and Goldstein [96] stated that near the step, the largest temperature 
difference was in the shear layer. As the reattachment point was 
reached, half of the temperature drop was across the shear layer and the 
other half was across the fluid near the wall. Sogin [98] found that 
for separation behind bluff bodies, most of the temperature drop was in 
the thin layer near the body. Vogel and Eaton [74] studied both the 
fluid dynamics and the heat transfer of flow through an asymmetric 
planar expansion. They found that the near wall region is important in 
determining the heat transfer rate. The sublayer dominates the heat 
transfer and is the reason the Stanton number and skin-friction 
coefficient are not well-correlated by the Reynolds analogy. This 
indicates that a constant turbulent Prandtl number is not correctly 
modeling the physical behavior in the recirculating region. They found 
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that the momentum thickness was much larger than the enthalpy thickness 
following reattachment. 

The reversed flow region affects the heat transfer well-beyond 
reattachment. The flow measured by Filetti and Kays [91] had not 
reached fully developed values afLer 14 step heights downstream of the 
step. Aung and Goldstein [96] said that their results indicated that 
the heat transfer values approached the flat plate values after 12 step 
heights. 

c. Numerical hydrodynamic Turbulent recirculating flows are 
much more diffi ult to predict than laminar recirculating flows. 
Turbulence models are well developed for attached boundary layers but 
experimental evidence has shown that the attached boundary layer 
assumptions are many times not valid in the recirculating flow. Most of 
the recent numerical predictions have used a variation of the k-e model 
of turbulence. Although relations for k and e can be rigorously derived 
from the Navier-Stokes equations, new unknowns are introduced that 
require modeling assumptions. These assumptions render the k-e 
transport equations approximations at best. Much of the computational 
research in turbulent recirculating flow has been carried out to improve 
the k-e model for recirculating flow. 

In the latter part of the 1970s, there was a sudden interest in 
numerical predictions of rapid expansion flow.* Briggs et al. [99] 
predicted the f low measured b’ Abbot and Kline [60] . Le Balleu- ^ud 
Mirande [100] and Kim et al . [101] u- 'd viscous -inviscid interaction. 


Ha Hinh and Chassaing [102J predicted flew through an asymmetrial 
expansion. Gosman et al. [103] predicted the flow through symmetric and 
axisymmetric expansions using primitive variables. Oliver [104] and 
Mehta [105] used vorticity-stream function variables. Taylor et al. 

[47] used a finite element method to solve the primitive variable form 
of the Navier-Stokes equations for flow over a backstep. He used a k- 
equation turbulence model with an empirical mixing length formula. 

Atkins et al. [106] predicted the flow through an asymmetric channel 
expansion. The above numerical predictions are discussed in more detail 
in Kvon and Pletcher [39]. 

Kwon and Pletcher [39] used viscous-inviscid interaction to predict 
the flow measured by Kim et al. [107]. They solved Laplace's equation 
in the inviscid region and the boundary layer equations in the viscous 
region. The FLARE approximation was used to march the boundary- layer 
equations through regions of reversed flow. They used the k-equation 
and an ODE for the length scale upstream of the step, and the k-equation 
and an algebraic equation for the length scale downstream of the step. 

Lokrou and Shen [108] solved a form of the boundary- layer equations 
by using a normalization of the velocity profiler to make them invariant 
in the streamwise direction. This reduced the PDFs to a system of ODEs. 
However, the theory fails near reattachment due tc flow curvature and 
instability. 

Sindir [109, 110] used four different turbulence models to predict 
the flow through asymmetric expansions with parallel walls and 
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nonparallel walls respectively. Two of the turbulence models were the 
k-E model and an algebraic stress model (which still requires solutions 
of the k and t equations). The other two were obtained by modifying the 
production term in the t -equation. He found the modified algebraic 
stress model superior in the reversed flow region. However, it 
predicted too slow a recovery after reattachment. He found that the 
best approach was to use the modified algebraic stress model in the 
reversed flow region and the nonmodified algebraic stress model after 
reattachment. The "best" model is thus regionally dependent. For the 
near wall region, he used the wall functions of Chieng and Launder 

[HI] - 

Hackman et al. [32] predicted the flow through an asymmetric 
expansion with two types of differencing schemes: (1) an upstream 

weighted difference scheme (UWDS) and (2) a skew hybrid upstream 
differencing scheme (SHUDS) for both Cartesian and curvilinear meshes. 
They solved the Navier-Stokes equations with the standard k-t turbulence 
model with law-of-the-wal 1 type wall functions near the wall. They 
found that the UWDS predicted shorter reattachment lengths, gave 
generally inferior predictions, and was much more sensitive to the mesh 
size. SHUDS was less grid dependent and gave better overall correlation 
with measurements. Hackman et al. thought that some of the poor 
predictions for turbulent flows was a result of the numerical scheme and 
not an inadequate turbulence model. Their computations over predicted 
the turbulent stress in the shear layer. This may be the reason for the 
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general under prediction of reattachment length by k-E turbulence 
models. Their model predicted a sharper return to low kinetic energy in 
the inviscid core than was measured. This effect is a common ailment of 
the k-E turbulence model for this type of flow. 

Srinivas and Fletcher [112] used a variation of an algebraic 
turbulence model to predict flow over the trailing edge of a flat plate 
and backward facing step. They used the finite element method to solve 
the compressible Navier-Stokes equations by a pseudo transient time 
marching technique. Since their wall pressure and maximum shear stress 
predictions were in good agreement with measurements, they argued that 
the algebraic eddy viscosity models predict most flow features well for 
wake and separated flow. 

Walterick et al. [73] predicted the flow over a backstep 
(asymmetric planar expansion) by solving the Navier-Si okes equations 
with a k-e turbulence model with the pressure fluctuation term of the k.- 
equation modeled in a nonstandard way. Their method predicted 
reattachment well. Their predictions using plug flow at the inlet gave 
shorter reattachment lengths than for the inlet condition with a 
boundary layer. 

d. Numerical heat transfer Chieng and Launder [111] used a 
modification of the TEACH-2E code to predict the flow through 
axisymmetric expansions. They tried the standard high-Reynolds -number 
k-E equation turbulence model with much attention given to a new set of 
wall functions, and a low-Reynolds -number k-E model. The low-Reynolds- 
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number model converged very slowly and predicted heat transfer rates 5 
times that of the measurements of Zemanick and Dougall [84]. The high- 
Reynolds -number turbulence model gave better predictions. However, in a 
later publication [113] it was reported that the original code had 
errors in it. When these errors were corrected, the predicted heat 
transfer was lowered by about 25%. The model was "fixed" by defining a 
variable turbulent Reynolds number related to the laminar sublayer 
thickness rather than assuming this turbulent Reynolds number to be a 
constant as in the first paper. This shows a common ingredient of 
turbulence modeling. Some models are so complex and contain so many 
constants and ad hoc functions, that these can be altered until the 
numerical predictions agree with experimental measurements. Even coding 
errors can apparently be overcome with appropriate turbulence modeling. 

Kang and Suzuki [114] computed the flow for a high speed jet in a 

pipe using the standard k-t turbulence model with law-of-the-wall type 

wall functions with constant values for c and Pr . They had to alter 

P t 

the enthalpy law of the wall to make their heat transfer rate 
predictions agree with experiments. 

Watkins and Gooray [115] and Gooray et al. [116] predicted the flow 
through asymmetric planar expansions and pipe expansions using a k-e 
turbulence model. The model was altered to include the effects of 



The 


44 


relations. They used a two pass procedure. The first pass was with a 
high-Reyno Ids -number k-E with wall functions to find the point of 
reattachment. The second pass consisted of the high-Reynolds -number k-e 
turbulence model upstream of reattachment and the low- Keynolds-number 
k-E equations of Jones and Launder [117] downstream of reattachment. 

The correlation between their predictions and experiments was very good. 
However, their expression for Pr^ for the Cartesian grid can be shown to 
be between 0.2 and 0.3 for fully developed equilibrium flow in a channel 
rather than the well accepted value of 0.9. 

Chieng [118] used a low-ReynolJs -number k-E turbulence model to 
predict the heat transfer in abrupt pipe expansions. Chieng' s 
predictions do not agree with the measurements of Zemanick and Dougall 
[84]. 

Amano [119] and Amano et al. [89] expanded on the two equation wall 
function method of Chieng and Launder [Hi]. Amano used a three zone 
wall function and did not require local equilibrium between production 
and dissipation in the E-equation. The predictions compare well with 
those of Zemanick and Dougall [84] for the d/D=0 . 54 expansion and high 
Reynolds numbei but do not compare as well for the d/D=0.43 expansion. 
The computations indicated that the maximum Stanton number was before 
reattachment and that the dependency of the level of heat transfer on Re 
is slightly less than that of the experimental data. T n Amano et al. 
[89], the computer predictions were compared with a concurrent set of 
experimental measurements. 


C. Scope and Contributions of the Present Study 


The purpose of this stud} is to determine the degree to which the 
boundary -layer equations can be used to model the flow through a region 
of flow reversal caused by an abrupt channel expansion. Since the 
previous studies that used the boundary- layer equations [39, 27, 42] 
provided only isolated comparisons with experiments and Navier-Stokes 
solutions, the limitations of the boundary- layer equations have not been 
previously defined. The purpose of the present work is to more clearly 
define the limitations of such solutions for both two-dimensional and 
axisymmetric expansion flows with respect to Reynolds number and 
expansion ratio. Determining the range of applicability of the 
boundary- layer equations is of practical importance since the effort 
required to solve the boundary- layer equations is an order of magnitude 
less than that required for the full Navier-Stokes equations. 
Furthermore, the constant property laminar boundary- layer equations are 
independent of Reynolds number. Therefore, the boundary- layer equations 
need to be solved only once for any given expansion ratio and the 
solution can then be applied through proper scaling for any channel 
Reynolds number. 

Figure 2 shows the geometry of the flows considered. Since both 
axisymmetric and symmetric planar expansions occur in applications, both 
types of geometry were considered. This allowed comparison with the 
results of as many other studies as possible. Only flows that were 
fully developed at the step were included in this work, so viscous- 
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Figure 2. Symmetric expansion geometry 

invisciu interaction was not applicable. Compressibility effects were 
negligible in the flow regimes considered so compressibility effects 
could be neglected. However, for _1 ows with heat trrnsfer, the 
variation of fluid properties with temperature was taken into account. 

Turbulence modeling of rapid expansion flows, especially with heat 
transfer, is a major challenge. None of the conventional algebraic or 
two-equation models work satisfactorily for the rapid expansion 
geometry. A number of model combinations were evaluated in this study. 
Turbulent flows were predicted by adding a turbulent viscosity to the 
molecular viscosity as was first proposed by Boussinesq [120]. The 
effect of various improvements and modifications of the k-e turbulence 
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model was considered. Some of the modifications tried included extra 
terms in the k and z equations, algebraic stress models, and variable 
turbulent Prandtl numbers. Special attention was given to different, 
methods of modeling the near-wall turbulence. The hydrodynamic and heat 
transfer predictions of the boundary- layer formulation were compared 
with other numerical results and with experimental measurements. 

The contributions of the present study are as follows. 

1. The limiting Reynolds number at which the laminar prediction 
begins to be poor is defined. 

2. The circumstances under which global iteration improves 
predictions are found. 

3. Those flow parameters well predicted by the boundary- layer 
formulation for laminar and turbulent flow are defined. 

4. Near-wall turbulence modeling was considered that did not 
assume a special form for velocity or temperature, i.e., the 
momentum, continuity, and energy equations were solved 
throughout rhe flow field including the near-wall region. 

5. A modification of the turbulent viscosity model of Johnson 
and King [121] was developed that gives improved results over 
the law-of-the-wall viscc ity model in the separated flow. 

6. This is the first study using the boundary-layer equations to 
predict turbulent flow with heat transfer through rapid 
expansions . 

7. An inexpensive numerical procedure is developed that has 
practical applications in predicting the flow and heat 
transfer in devices having large regions of separated flow. 
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II. GOVERNING EQUATIONS 

In this chapter, the theoretical background needed to solve the 
problem of flow and heat transfer in a sudden expansion is developed. 

The relevant partial-differential equations will be developed with their 
associated boundary and initial conditions. The compressible, variable 
property, turbulent boundary -layer equations will be developed first. 
These equations are the most general needed in this thesis. Particular 
simplifications of the compressible turbulent boundary- layer equations 
will be made as needed. It is a simple matter to make the necessary 
simplifications of the general turbulent boundary- layer equations to 
obtain, for instance, the laminar incompressible equations. A special 
form of the laminar incompressible equations can then developed that is 
independent of the Reynolds number. The turbulence models used will 
then be discussed. Finally, the relevant engineering parameters are 
presented . 

A. Variable Property Turbulent Boundary-Layer Equations 

The continuity equation is by far the easiest equation to derive. 
For this reason, this equation will be developed in more detail than the 
momentum or energy equations in ordei to show the methodology and 
introduce the necessary background information. 
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1 . Conservation of pass 

A mass balance of the i-'iid passing through a fixed control volume 
produces the following equation 

|f + V*( P V) = 0 (2.1) 

Equation (2.1) is valid for any nonreacting continuum of fluid, 
including fluid with turbulence. However, to discretize the equation 
fjr computer simulation of the turbulent motion, a grid fine enough to 
capture the very small length scales of the turbulent motion must be 
>.sed. For a practical problem, this would require the solution of a 
number of algebraic equations that is beyond the capabilities of a 
present day computer. Since turbulent motion is characterized as random 
motion in which statistical averages are meaningful [122], time averaged 
equations can be used in nlace of the instantaneous set. To average the 
equations, the instantaneous velocity of the fluid certain 

properties) is considered as the sum of an averaged value and its 
fluctuating value [123]. 

f = f + f' (2.2) 

Hare, f represents a velocity component or property of the fluid and is 
defined as 

f =» / t+At f(t)dt 
t 


(2.3) 
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Note that by definition f' = 0.0, and f *g =0.0 where g is some averaged 
variable. The time of averaging At is long enough to average out the 
fluctuations due to turbulence but not long enough to influence the 


variation of the mean flow with time. The density and velocity vector 
are written in the following form 

P = P + p' 

(2.4) 

V = V + V’ 

(2.5) 

Equations (2.4) and (2.5) are substituted into (2.1) and 

the equation is 

time averaged. After canceling the terms that are zero, 
equation results: 

a time averaged 

V • (pV + J/V* ) = 0 

(2.6) 

It is now convenient to introduce the following type of 
is often used for gas mixtures [123] 

averaging that 

: = f + f" 

where 

(2.7) 

f = 

p 

It can be shown that 

(2.8) 

» 1 _ — —If 

p u = - pu 

(2.9) 
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Using (2.9) in (2.6) and making the assumption of steady flow, Eq. (2.6) 
becomes 

7 • (pV) = 0 (2.10) 

Under the boundary-assumptions [123], 

u = u (2.11) 


Equation (2.10) is now written in expanded form for two dimensional 
planar and axisymmetric flow with u substituted for u. 

|^(r m pu) + fy(r m pv) = 0 (2.12) 

Equation (2.12) is valid for compressible flow for both two-dimensional 
Cartesian and axisymmetric geometries. When the geometry is Cartesian, 
m is taken as zero; when axisymmetric, m is taken as one. 

It is convenient (but not essential [39]) when using the present 
scheme to employ the strean f motion, , in place of v as a variable. 
This is done merely to make the equations more suitable for external 
flows with separation. ."hen using the boundary- layer equations for 
external flows with separation, the displacement thickness is specified 
as a boundary condition to overcome the singularity at separation. If 
u-v variables are used, the displacement thickness must be obtained by 
integration of the velocity profile; if u-V> variables are used it is 
obtained directly without integration. Au expression for the stream 
function is 


M - . 


r pv 


(2.13) 
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Substituting Eq. (2.13) into the continuity Eq. (2.12) gives a new form 
of the continuity equation 


M = 

3y 


m 

r pu 


(2.14) 


Note that the definition of i> for compressible flow corresponds to a 
mass flow rate rather than a volumetric flow rate as for incompressible 
flows . 

Equation (2.12) is the same as the continuity equation for laminar 
flow neglecting the tilda and bar notation above the averaged variables. 


2. Conservation of momentum 

The momentum equation is merely a statement of Newton's second law 
for fluid r lowing through a nonacceierating control volume. In vector 
form for compressible flows, it is [124] 

= - Vp - [V*-T] + pg (2.15) 


Note that T is the stress tensor and g is the gravity vector. To obtain 
the equation valid for compressible turbulent boundary layers, a 
derivation similar to that followed for the continuity equation must be 
followed. However, for the momentum and energy equations this is so 
lengthy that the steps will only be outlined here. (See Cebeci and 
Smith [123] and Anderson et al. [4] for the details of the derivation.) 

The main steps required to obtain the compressible turbulent 
be undary- layer momentum equations are as follows. 
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(1) Substitute the sum of the averaged and fluctuating value for each of 
the instantaneous velocity components, density, and pressure. 

(2) Time average the equations and cancel the terms that are zero due to 
the time averaging. 

(3) Substitute the continuity equation into the momentum equation. 

(4) Neglect the body forces and the mean or averaged unsteady terms. 

(5) Make the boundary- layer assumptions and cancel the higher order 
terms by doing an order of magnitude analysis. The main assumptions 
made are: 


3 3 _ 

— »— ,u»v,u = u 
3y 3x ’ 


For subsonic flow, the x-momentum boundary- layer equation is then 


— 3u , 3u dp , 1 3 m, 3u 

pu si + pv 3y = ix * f Ty r (v Vy ' 

The y -momentum equation reduces to 


dp 3 r- II II. 

' 5? ' ^ (pv v > = 0 


(2.16) 


(2.17) 


The pressure gradient in the y-direction is assumed to be small compared 
to the one in the x-direction, so the y-momentum equation can be 
neglected [4J. 

Equation (2.16), neglecting the bars indicating a time averaged 
quantity and the Reynolds stress term, pu'v' , is the same as the x- 
momentum equation governing laminar boundary- layer flow. The extra term 
(-pu'v') arises due to the convection of turbulence. It is grouped with 
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the laminar-like viscous term because u 1 v ' can be thought of as 
increasing the stress through exchange of momentum in the fluid. Due to 
the -pu'v' term, there are presently more unknowns than equations. The 
extra equations needed to solve the syste 1 are obtained through 
turbulence modeling. 

At this point, it is convenient to introduce a turbulent viscosity 
defined as 

\ = - P uV < 2 - 18 > 

as vas first suggested by Boussinesq in 1877 [122]. The Boussinesq 
assumption is actually the first turbulent modeling assumption made. 
Further turbulence modeling assumptions and the evaluation of the 
turbulent viscosity is discussed in Section C of this chapter. Of 
course if the flow being predicted is laminar, all turbulent fluctuation 
velocities are zero, so is also zero. 

At this point, it is convenient to drop the bars and tildas in Eqs. 
(2.12) and (2.16). By substituting Eq. (2.18) into Eq. (2.16), the 
momentum equation becomes 


3u . 3u _ dp 1 3 m, . .3u 

pU to + to + ^ 3^ 1 


(2.19) 


The unbarred variables are recognized to represent time averaged 
quantities in turbulent flow. 


Equation (2.19) is a parabolic equation in x. It is normally 
solved by marching from a starting position to the desired location in 
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the x-direction. As the solution is marched down the channel for u > 0, 
the solution should only be influenced by the domain behind the x- 
position that has been reached by marching. What is beyond the point 
reached by marching in the x-direction ^downstream) should not influencT 
the solution of a truly parabolic equation, or the formulation becomes 
unstable. However, with flow reversal, flow moving in the negative x- 
direction will influence the solution from downstream. This downstream 
influence causes the streamwise marching solution of the parabolic 
equations to be unstable unless special measures are taken. The problem 
lies in the x-convective term u3u/3x since it is the term that 
contributes the downstream influence. 

Reyhner and Flugge-Lotz [53] proposed that the convective term in 
the x-direction be replaced by c|u|3u/3x where c is a small positive 
constant or zero. This approximation is referred to as the FLARE 
approximation. For rapid expansions, the velocity in the reversed flow 
is about 10% to 20% that of the main flow stream, so this assumption 
seems valid. For the momentum and energy equations, c will be taken as 
one if u is positive; c will be taken as zero if u is negative. Thus, 
any downstream influence to the parabolic equations is cancelled and the 
formulation is stable. Governing partial-differential equations for 
variables associated with turbulence will be developed in following 
sections that have an x-convective term similar to Eq. (2.19). In 
solving these turbulence modeling equations, taking c as 0.0 when u was 
negative caused numerical instabilities near the wall associated with 
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round-off errors in the computer. For stability reasons, c was set to 
1.0 in the turbulence modeling equations when using the FLARE 
approximation . 

Repetitive global iteration through the flow field can be used in 
place of the FLARE approximation to provide stability to the boundary- 
layer equation marching procedure. The x-convective term can be 
approximated with a local-upwind finite-difference as opposed to 
disregarding or approximating this term as is done with the FLARE 
approximation. 

We can also use the previously introduced stream function to 
eliminate v from Eq. (2.19). After making the FLARE approximation and 
using Eq. (2.13) to eliminate v from Eq. (2.32a), the resulting x- 
momentum equation valid for turbulent compressible flow is 


cpu 


9u 1 3^ 3u dp 1 3 m. .3u 


9x m 9x 9y dx m 9y 
r J r J 


t 9y 


( 2 . 20 ) 


Equation (2.20) is singular at the centerline for axisymmetric 
geometries. When r = 0, the y-convective term can be removed from Eq. 
(2.20) since v is zero and no convection in the y direction will occur. 
l'Hospital's rule can then be used to find the valid representation of 
the diffusion term as r approaches zero. The momentum equation valid 
for axisymmetric flow with r = 0 is 


cpu 


9u _ _dp 
9x dx 



<“+v| 


( 2 . 21 ) 
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3. Conserv it ion of energy 

The energy equation is simply a statement of the first law of 
thermodynamics for fluid passing through a stationary control volume. 
For flows */hich are not isothermal, the energy equation must be solved 
along with the continuity and momentum equations. There are two main 
reasons why the energy equation must be introduced: (1) heat transfer 
quantities or temperature fields are desired or (2) the change of the 
fluid properties with temperature will affect the hydrodynamic 
predictions . 

The energy equation can be written in terms of the total enthalpy, 
H. 


|^CpH) + V(pHV) = |2 + V{(t*V] - q} 
K is defined as 


( 2 . 22 ) 


H = h + 2U.U. 


(2.23) 


The stress tensor, t, depends on the coordinate system (see Anderson et 
al. [4] for the different forms). The heat flux vector, q, can be 
written in indicial form as 



The body forces have be *n neglected. 

Cebeci and Smith 1123] show the details of obtaining the form of 
the energy equation used in this thesis. An outline of the steps 
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required to obtain the desired form of the energy equation follows. 
Equation (2.22) is time averaged in the same way as the continuity and 
momentum equations discussed in the previous sections. The boundary- 
layer assumptions are then used to eliminate the higher order terms 
through an order of magnitude analysis. Under boundary- layer 
assumptions, H = H. For steady, subsonic, compressible flow, the energy 
equation becomes 


3H . 3H 
pu— + pvr— 
K 3x K 3y 


1_ 31 m/y_ 3H 
^m 3yl r \Pr 3y 


pc v'T’ + uf(l-~)p~-pu'v 


r Pr 3y 


(2.25) 


Equation (2.25) is the same as the laminar energy equation except 
for the two extra terms involving the fluctuating turbulent quantities, 
-pCpT'v' and -pu'v'. The heat flux caused by the T'v' term is assumed 
to be proportional to 3T/3y. A turbulent Prandtl number, Pr , is 
defined by the following equation 

(2 - 26) 

Pr^ must be determined by turbulence modeling. When Eqs . (2.18), 

(2.26), and (2.13) are. substituted into Eq. (2.25) and after making the 
FLARE approximation in the same way as was done for the momentum 
equation, the turbulent compressible boundary- layer energy equation 
becomes 


3H 

1_ 3£ 

IS 


3x 

m 3x 

3y 



r 


1_ 

[!_ "/I 

v_ + y t 1 

15 1 

mi 

r 

lay \l 

Fr PrJ 

5y 1 



+ P. 



(2.27) 
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The bars have been dropped from Eq. (2.27) and all the flow variables 
are understood to be averaged values. 

Equation (2.27) is singular when r = 0 for axisymmetric geometries. 
By removing the y-convective term and using l'Hospital's rule, the 
singularity can be removed. Equation (2.27) becomes 


3H 

pU I^ = 


2 *#- + 
Wr 




yd-fe) 


+ y* 



(2.28) 


Equation (2.28) is the governing momentum equation when r = 0 for 
axisymmetric geometries. 


4. Boundary and initial conditions 

a. Boundary conditions Two boundary conditions for both u and 
H were used, one at y = 0 and one at the centerline. One boundary 
condition at the wall was specified for 4 >. In addition to the wall 
boundary condition for another restriction the solution must satisfy 
is the channel mass flow constraint. The mass flow constraint is merely 
a mass balance across the entire channel cross section. The channel 
mass flow constraint is not actually used as a boundary condition but 
will be necessary later in order to solve for dp/dx. 

The no-slip condition at the wall and the symmetry condition at the 
channel centerline give the boundary conditions for the velocity as 


dll 

u(x,0) = 0 , — (x,y . . ) = 0 

8y centerline 


(2.29) 


The stream function boundary condition and mass flow constraint for 


the channel cross section are respectively 



♦(*•0) ■ 0 ■ *t*-Werline ) = ' D/ V m u( 0 , y )dy (2.30) 

The boundary condition for H at the centerline is the symmetry 
condition 


3^ (x * y 


centerline 


) = 0 


(2.31) 


Two different H boundary conditions are possible at the wall depending 
on whether the temperature or the heat transfer r^e is specified. For 
specified wall temperature, T^(x) , the boundary condition is 


H(x,0) = c T (x) (2.3a) 

p w 

For specified heat flux, q^x), the boundary condition is 

-kfjdf-f > " <yi*Cx) (2.32b) 

b. Initial conditions The initial conditions depend on whether 
the flow is turbulent or laminar. Above the lip of the step (y > h) , a 
fully developed turbulent or fully developed laminar inlet profile was 
used as the inlet condition for u. Along the face of the step, the 
logical condition due to the no slip requirement is 


u(0,y) =0, 0 1 y < h 


(2.33) 


Howe' Acrivos and Schrader [27] argue that this is not correct for 
the boundary -layer equation set. They say that a nonzero velocity 
should be used to take into account the effect of the fluid returning 
from downstream in the recirculating region. Acrivos and Schrader used 
the following initial condition on the face of the step. 
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u(0,y) = u c (y), 0 S yS h/2 (2.34) 

u(0,y) = - u c (h-y), h/2 < y <: h 

The fluid velocity was taken as zero for this study except for in a few 

cases used to determine th« ef 4 et of the nonzero velocity on the face 
of the step. 

Once the initial velocity is decided upon, the following integral 
gives 4* 

^(0,y) - / y pr m u(0,£)d£ (2.35) 

0 

where ( is a dummy variable of integration. 

Since the experimental studies compared with used an insulated step 
face and a fully developed temperature profile, the inlet condition for 
H is 

f~(0,y) =0, 0 < y < h (2.36a) 

H(0,y)= c p T(0,y) + ^(O.y), h < y < D/2 (2.36b) 

For a fully developed inlet profile and an unheated entry length, 

T(0,y) = T b (0) (2.3V) 

Where T, (x) is the mean or bulk temperature. 

D 

5. Equation of state 

In order to solve the above conservat-' on equations, p, y, k, and Pr 
must be specified. This is done by introducing an equation of • ate. 
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In general, any property of the fluid, f, is a function of temperature 
and pressure, i.e. , 


f = f(T,P) 


(2.38) 


Since one of the main purposes of this study is to compare with 
experimental measurements, the fluids of interest are water and air at 
moderate temperatures and nearly atmospheric pressures. For these 
restrictions, all the properties of water are very weak functions of 
pressure (125]. The density of water for all practical purposes is 
constant. For air near atmospheric pressure, u, k, and Pr are all very 
weak functions of pressure. For air, the density is found from the 
ideal gas equation 


= E_ 
RT 


(2.39) 


where R is the gas constant. Thus, all properties except the density 
for air vary only with temperature and can be expressed as 


f * f(T) 


(2.40) 


The particular functions used to approximate the fluid properties as a 
function of temperature for air and water are given in Appendix A. 


B. Laminar Constant Property Nondimensionalized Equations 


For laminar constant property flows, the boundary- layer equations 
can be nondimens ionalized such that they are independent of the Reynolds 
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number. The method of nondimensionalizing has been previously used for 
developing flow in channels, it reduces the hydrodynamic solution to a 
function of the size of the expansion, the inlet conditions, and the 
boundary conditions. 

For constant density, p can be removed from the partial- 
differential operators of Eq. (2.12) and canceled from the equation. A 
different stream function must be defined for the incompressible case 
than was used for the compressible case in the previous section. The 
stream function for incompressible flow is defined as 


v = * 

3x cp 


The constant property continuity equation can be written as 


(2.41) 


U = 3? ( V 


(2.42) 


To '.ondimensionalize the conservation equations, the following 
dimensionless variaMu are introduced: 


U = - , X = , Y * * n = 

u. dRe d ’ u* 

l i 

jb 

. _ cp _ r - _ dRe dp 

” » * j * P ~ H 


,1-Hn 
d u 


pu . dx 


(2.43a) 

(2.43b) 


Again, m is zerc if the geometry is two-dimensional; m is one if the 
geometry is axisymmetric . The Reynolds number is based on u^ an 1 d. 

Using Eq. (2.43) in Eq. (2.42) gives the dimensionless continuity 
equation 


U , 

R" 3Y 


(2.4 
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Substituting Eq. (2.43) into tue laminar form of Eq. (2.16) gives the 
aimens ionless x-momentum equation 



i_ II _ o . i_ L. rR m au 
R m 3X 3Y P R m 3Y QK 3Y 


(2.45) 


Equations (2.43) with the laminar energy equation (Eq. 2.27) gives tne 
dimensionless energy equation 



I_ II in = 1. 3_/ R ail_ + n 

R m 3X 3Y R m 3Y\ K iPr 3Y u 



Note that the FLARE approximation has been made i. 
for the dimensional equations in the previous sect?' .. 
The bour.Hary conditions ::or U arc 


(2.46) 
way as done 


U(X,0) = 0 , |^(X, 1 + |) = 0 


(2.47) 


The stream function is set to zero at the caannel wall (Y(X,0) = 0). 

The stream function at the centerline is determined from the inlet 
Dro. at the step. For a parabolic, fully developed inlet 

Y(X,1 + 4) = planar expansion (2.48a) 

Y(X,1 + ^) = |, axisymmetric expansion (2.48b) 

The boundary condition for ti at the centerline is the symmetry 
condition 


|i(X,l + J) = 0 


(2.49) 
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There are two possible boundary conditions for ti at the wall depending 
on if the heat flux or the wall temperature is specified. If the wall 
temperature is specified, the condition is 


U(X,0) = Jyc T (X,0) 
i P 

If the heat flux at the wall is specified, the condition is 


-|f(x,0) - ♦ If jfl 


y=0 


where is a dimensionless wall heat flux given by 
Pe % 


puT 


(2.50a) 


(2.50b) 


(2.51) 


Pe is the Peclet number (?e = Re * Pr) . 

A parabolic initial condition for U at the inlet uas specif.; vi. U 
was taken as i’ero along the face of the step. The initial condition for 
the stream function is obtained by integrating U from the wall to the 
centerline as follows 


T( 0 ,Y) = /Vu(0,$)dC 

0 

The initial condition for the nondimens ional total enthalpy is 
c T(0,Y) 


(2.52) 


H(0,Y) = -2—r 


+ |u*(0,Y) 


(2.53) 
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C. Turbulence Modeling 

In a previous section, Eq. (2.18) defined a turbulent viscosity, y t 
according to the Boussinesq approximation. Equation (2.26) defined a 
turbulent Prandtl number, Pr . In order to predict turbulent flow, y 
and Pr^ must be approximated by turbulence modeling. The modeling used 
in this work is discussed in this section. 

2- Equilibrium tu rbulence equations 

Turbulent boundary -layer flow with no separation can be divided 
into two main regions: (1) the inner region which is not highly 

influenced by the pressure gradient, and (2) the outer region which is 
highly affected by the pressure gradient. The inr~r region can be 
further divided into three parts: a laminar sublayer, a buffer region, 
and a fully turbulent region. In the laminar sublayer, the molecular 
viscosity dominates; in the fully turbulent region, the Reynolds 
stresses dominate. Both the stresses due to molecular viscosity and 
Reynolds stresses are important in the buffer region. 

Dimensional analysis shows that a nondimensioral velocity, u + , and 
a nondimensional distance from the wall, y + , are important in describing 
the flow in the inner region. They are defined as 

+ _ U t yp + _ u (2.54) 

y V ’ " u t 

1/2 

where u^ is a turbulent velocity scale (t^/p) and * * s the shear 
stress at the wall. When y + and u + are used to plot turbulent boundary- 
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layer velocity profiles with no separation, the profiles collapse into 
one curve in the inner region. Experiments have shown that the time 
averaged thicknesses of these three layers for smooth walls ..re 
approximately 

laminar sublayer 0 < y + < 5 

buffer layer 5 < y + < 30 

fully turbulent layer 30 < y + 

The molecular viscosity of gases can be calculated by 

= p( length scale) (velocity scale) (2.55) 

Prandtl adopted this idea for turbulent ilow and proposed a length 
scale. It, for the inner region that is proportional to the distance from 
the wall. Van Driest later modified l by multiplying it by an 
experimentally determined damping function D [126J. The modified length 
scale becomes 


It = xDy 


(2.56) 


where ^ = 0.41. The van Driest damping function, D, is given by 


D = 1 - exp(-y + /A + ) 


(2.57) 


where A T is usually taken as 25 or 26. Prandtl proposed a velocity 
scale as 


(velocity scale) 



(2.58) 


The turbulent viscosity for the inner region can then be expressed as 
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,2 |3u 

w t = ** m 


(2.59) 


Equation (2.59) is known as the Prandtl mixii g length model for 
turbulent viscosity. With specified, che governing equations can be 
integrated for the inner region. 

By making the Couette flow assumption (3u/3x is very small and can 
be neglected) the momentum equation is reduced to an ODE [127]. This 
assumption is valid near the wall. By neglecting the turbulent 
viscosity, this ODE can be integrated in the laminar sublayer cf the 
inner region to show that 


+ + 
u = y 


(2.60) 


This is equivalent to saying that the shear stress is constant in the 
near-wall region. Similarly, by neglecting the laminar viscosity in the 
fully turbulent part of the inner region and using Eq. (2.59), this ODE 
can be integrated to show that 


u + = ~ln(y + ) + B 


(2.61) 


where B is a constant near 5.0. The region for which Eq. (2.61) is 

valid is some es referred to as the logarithmic layer. Equation 

(2.61) is often called the "law of the wall". The velocity in the 

buffer region, u + = f(y + ) must be obtained experimentally. 

There have been attempts to modify Eq. (2.56) to make the mixing 

length more applicable for reversed flow. Reeves [123] and 

McD Galbraith and Head [129] recommend that Eq. (2.56) be multiplied by 
1/2 

(t /x ) The n-lxing length for the inner region then becomes 

max w 
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x 

t = k Dy (2.62) 

w 

Carter and Wornom [130] and Pletcher [131] recommend that the van Driest 
damping function, D, by based on the maximum shear stress in the 
boundary layer instead of the wall shear stress. The modified van 
Driest damping factor is then 


D = 1 - exp[-( f |!^| ) 1/2 y/A + ^ 

yl 3ylmax 3/ 

For an attached boundary layer, Eqs . (2.63) and (2.63) reduce to Eqs 


(2.63) 


(2.56) and (2.57). 

The value of 25 for A + corresponds to a boundary layer with a zero 
pressure gradient. Kays and Moffat [132] proposed an empirically ~ased 
function for A + valid for nonzero pressure gradients as 


A + = 25 . 0/ (ap + + 1) (2.64) 

where p + is given by 

. V (dp/dx) 
p = 1/2 3/2 

p T 

W 

and 

a = 30.2, p f < C 
a = ,?0.6, p + > 0 

Johnson and King [121] expressed the turbulent viscosity according 

to Eq. (2.55) but used the square root of the maximum Reynolds stress, 

— — 1/2 2 

(u'v' ) ' , is the turbulent velocity scale. They used D 1c Ky as the 
m 15 

length scale where is the van Driest damping function based cn an A 
value of 15 instead of the usual 25. The expression for V* t then becomes 


,'W 
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W t = P D l5 Ky( " U,V, m )1/2 (2-65) 

Equation (2.65) can be used in place of Eq. (2.59) for the turbulent - - 
viscosity expression in the inner region. 

For the outer region, £ is often taken as proportional to the 
boundary layer thickness with Eq. (2.58) used as the v - ocity scale [4]. 
Another common method is to use a constant times the displacement 
thickness as the length scale and the velocity at the outer edge of the 
boundary layer as the velocity scale. The expression for the turbulent 
viscosity in the outer region then becomes 

P t = p(K6*)(u e ) (2.66) 

* 

where K is a constant (- 0.016), 0 is the displacement thickness as 
evaluated for incompressible flow (145), and u^ is the velocity at the 
edge of the boundary layer. 

Equations (2.59), (2.65), and (2.66) describe zero equation 
turbulence models. One way of classifying turbulence models is to add 
the number of PDEs used in the turbulence model. Each PDE counts as one 
equation; each ODE counts as one-half an equation; each algebraic 
equation counts as zero equations. 

Zero equation turbulence models are not able to take into account 
the effects of diffusion or convection of turbulence length scales or 
velocity scales. They only balance the production and dissipation of 
turbulence quantities [3] . To account for convection and diffusion, 

PDEs are required to introduce directional rates of change. When the 


71 


dissipation of turbulence balances the production, the turbulence is 
said to be in local equilibrium. Therefore, zero or algebraic 
turbulence models are valid for regions in local equilibrium. In the 
shear layer for flow over a step, turbulent kinetic energy production 
exceeds dissipation [24]. This turbulent kinetic energy must, be 
convectfed and diffused away. This indicates that in the region affected 
by the shear layer, a zero equation model is inadequate a-d models 
involving one or more PDEs or ODEs must be used. 


2 . k-c equations 

In this section, the two-equation turbulence models investigated in 

this work are described. For the k-e model, turbulent length and 

velocity scales for use in Eq. (2.55) are described by PDEs [3]. The 

1/2 

velocity scale used is k , where k is called the turbulent kinetic 
energy and is defined as 

k = ~ (u’ 2 + v’ 2 + w’ 2 ) (2.67) 


The length scale is defined as 



( 2 . 68 ) 


wh°re Cp is a constant listed in Table 2 and t is the dissipation rate 
of k. The expression for the turbulent viscosity is given by 


, . __ is. 

p. = pc - 
t p e 


where c is usually 0.09 
V 


(2.69) 
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Table 2. Turbulent Modeling Constants 


Constant Value Reference Constant Value Reference 


c 

U 

0.09 

[133] 

a 

2.2 

[134] 

C D 

0.164 

[39] 

y 

0.55 

[134] 

C 1 

1.35 

[135] 

C 01 

3.0 

[115] 

C 2 

1.8 

[135] 

C 82 

0.33 

[115] 

C 3 

0.0115 

[135] 

c eiw 

0.5 

[115] 

C 4 

0.5 

[135] 

c 

w 

2.44 

[115] 


At this poii. . , the unknown has been expressed in terms of two 
other unknowns, k and z. Transport PDEs can be derived for both k and z 
from the Navier-Stokes equations, but these PDEs involve other unknowns. 
The formulation must be "closed" by modeling assumptions to provide the 
necessary equations for the additional unknowns . After these modeling 
assumptions, the k and e equations can be written as 


-Mf* f* - i Jjl*r $1 ♦ s 


(2.70) 


where 0 is either k or z , depending on the equation desired. 3, is a 

9 

source term. 

Two different sets of source ter- and T's were used in this study: 


a high Reynolds number form used by Launder and Spalding [133], end a 
low Reynolds number form developed by Chien [135]. Chien added extra 
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terms to T and S to make the original higl. turbulent Reynolds number 

0 v 

equations used by Launder and Spalding applicable in regions of low 
turbulent Reynolds number flow near the wall. 

In referring to turbulence models, the Reynolds number of 
turbulence or turbulent Reynolds number is not the same as the Reynold 
number based on the the channel width or pipe diameter and the average 
inlet velocity. The dividing line between the low and high turbulent 
Reynolds number regions is not universally or clearly defined for all 
types of flow. There are several ways to define the turbulent Reynolds 
number depending on the velocity scale and the length scale used. If y 
is the length scale and u^ the velocity scale, the turbulent Reynolds 
number is merely y + . For this definition, the high turbulent Reynolds 

number flow for an attached boundary layer would be that part of the 

flow for which y + > 30. The low Reynolds number flow would correspond 
to y + < 30. Hereafter, the high and low turbulent Reynolds number k-e 
models will be referred to as the nigh-Reynolds -number k*e model and the 
low-Reynolds -number k-e model respectively, realizing that the Reynolds 
number referred to in these cases is the Reynolds number of turbulence. 

The expressions for V, and S. for the k and e equations for both 

0 0 

the high and low Reynolds number cases are given in Table 3. Those 

terms within the dotted vertical lines are those added by Chian in his 

low Reynolds number model. Chien -llso found it necessary to modify the 
turbulent viscosity relation by multiplying Eq. (2.69) by an empirical 
function as follows 

k 2 + 

\ = I (1 '“ p( * c 3 y )] 


(1 ) 


Table 3. 

Expressions 

for S, and V. for the k-E turbulence model 
0 0 

0 

r 0 

S 0 

k 

p + P t 

*»<$*>■ 

£ 


C -u f — ) 2 - p^fc 'f'e +'^ e ( " c 4 y+)| l 
1 rt 3y p k l °2 ' « e 'py* |J 



(f = 1 - 0.222e^ (’P k */6pe) a ]j 


The boundary conditions for the low-Reynolds -number k-e equation 
are applied at the wall as 

k = e = 0 (2.72) 

The high-Reynolds-number k-e equations a^e not valid near the wall. 
Boundary conditions must be applied away from the wall or wall functions 
must be used to approximate k and e near the wall [4]. 

Two expressions for k and £ in '.he near-wall region were used in 
this study when solving the high Reynolds number k-£ equations. One is 
based on the turbulent viscosity near the wall as expressed by Eq. 

(2.5?' £.nu the other is based on as given by Eq. (2.65). 

The wall model based on Prandtl's mixing length model was derived 
as fellers. For flow in equilibrium (production = dissipation), the 
convective and diffusive terms of the high Reynolds number case of Eq. 
(2.71) are neglected. Inserting the expression for the mixing length 


gives 
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e 2 iv 

May) “ pC D L 


3/2 


(2.73) 


For this equilibrium case and the constants given in Table 2, it can be 
shown that [3] 


2 


Using Eq. 


2 


2 


e 3/2 /c n L = 1.002 L « L 
y ' D 

(2.59) to eliminate y t gives 

3u| 3 _ k 3/2 

3yl C D i 


Solving for k gives 


k = c' 2 / 3 l 2 A 2 

D C 3y' 


(2.74) 


Equation (2.74) can be used to find k in the near-wall region and as a 
boundary condition for the k-equation PDE. Once k is known, e is found 
from 

3/2 


£ = c d'l 


(2.75) 


Equations (2.74) and (2.75) should be used in place of the PDEs 
expressing k aiid e for y + < 30 to 100. 

Simpson et al. [76] argues that measurements in a separation bubble 
on a flat pJate caused by pressure gradients indicates that that che 
law-of-the-wall model (which is based on Eq . 2.54) is not valid even at 
the wall in the separated flow. He states that the flow is dominated by 
turbulent fluctuations that are comparable in magnitude to the mean 
velocities. However, several investigators have used the equilibrium 
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law-of-the-wall model for the near-wall flow in conjunction with two- 
equation models away from the wall with good success [32, 116]. 

The near-wall expressions for k and t based upon the inner 
viscosity model of Johnson and King [121] were derived by first 
substituting Eq. (2.65) into the equilibrium expression given by Eq. 
(2.73). After assuming L is given by KDy and solving for k gives 


k ■ g>V»v >) l ' 2 i 2 ' 3 

Once k is known, the expressions for given by Eqs. 
are equated and solved for k to give 

, 2 . rn 2 , 7 r -V 1 / 2 ^ 

e = c/ /[D l 5 kyC-u v _) ] 


(2.76) 

(2.65) and (2.69) 


(2.77) 


Equations (2.76) and (2.77) provide alternative expressions for k and e 
in the near-wall region. 


3. Algebraic stress model 

The algebraic stress model (ASM) is becoming a popular variation of 
the k-e turbulence model. In the ASM, an expression for c^ is derived 
that is used in place of the constant value of 0.09. The functional c 
is used in conjunction with the k-E equations and Eq. (2.69) to 
calculate the turbulent viscosity. The ASM includes some effects of 
pressure-strain interaction and streamline curvature. 

The algebraic stress model used is that developed by Rodi [3]. 

Exact PDEs can be derived for the six different Reynolds stresses of 
which only u'v' is significant for boundary- layer flow. After modeling. 
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these transport equations for the Reynolds stresses contain gradients of 
the Reynolds stresses only in the convection and diffusion terms. Rodi 
assumes that the transport of ujuj Is proportional to the transport of 
k. An equation for the transport of k can be obtained from Eq. (2.70) 
as 


ft - Diff(k > ■ \ 


where D( )/Dt represents the toctl derivative, and Diff( ) represents 
the diffusion operator. Rodi assumes the proportionality constant 


between the transport of u!u! and the transport of k as u^u-/k. The 

^ J 1 

following equation can then be written for the transport of the Reynolds 


stresses 


- ST^/k t^-Diff(k)] 


(2.78) 


Equation (2.78) is used in the Reynolds stress transport equations with 
modeling assumptions to derive algebraic equations for the Reynolds 
stresses. Thxs is possible since the gradient terms have been removed 
by the proportionality assumption expressed as Eq. (2.78). For thin 
shear layers with no bouyancy terms included, u'v' is 


-u'v' = c 


k 2 3u 
V £ 3y 


(2.79) 


(an eddy viscosity relationship) where c 


V 


is now a function given by 


_ 2 n y , a-l+YP/t 

c u 3 (1 ' y) 2 
V J (I-l+P/e) 


(2.80) 


and P is 
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D l,3u,2 

P ■ "t p<5P> 


The values of a and t are given in Table 2, 


4. Turbulent prandt 1 number 

When solving the energy equation for turbulent flow, the effects of 
turbulent mixing are included in the turbulent Prandt 1 number, Pr fc , as 
defined by Eq. (2.26). The Reynolds analogy [127} states that the 
turbulent diffusivity for momentum, p^/p, is equal to the turbulent 
thermal diffusivity, k^/Cpc^), where k t is an effective thermal 
conductivity coefficient. This is based on the idea that if the 
turbulent motion of the fluid is dominant over molecular diffusion, the 
heat and the momentum turbulent diffusion should proceed at nearly the 
same rate since the same mechanism is responsible for both. For 
molecular diffusion, the Prandtl number is equal to the molecular 
diffusivity of momentum (kinematic viscosity) divided by the molecular 
thermal diffusivity. If the turbulent Prandtl number is describe'* in 
the same way, it would equal one, since the two dif fusivities are equal 
under the Reynolds analogy. 

Experimental measurements have shown that for molecular Prandtl 
numbers greater than 0.5, Pr is closer to 0.9 if it must be a constant. 
However, Pr c in be as high as 2.0 near the wall and as low as 0 8 in 
the fully turbulent region [127]. Pr^ = 0.9 seems to be an average 
value that works well for aii . Several recent turbulent heat transfer 
predictions have used P .• = 0.9 with reasonable success [111, 115]. 
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It is possible to derive expressions for Pr^ in a method similar to 
that used to derive an expression for c^. First, transport equations 
for the three turbulent scalar fluxes u^T 7 are derived. The gradient 
terms in these equations are removed by modeling assumptions to leave an 
algebraic expression for Pr t> Gibson and Launder [136] assumed that the 
transport of the scalar fluxes was proportional to the transport of k as 
Rodi did to derive an expression for c^. 

Watkins and Gooray [115] derived an expression for Pr^ in 
streamline and Cartesian coordinates. Their derivation was for a 
general elliptic problem. When the boundary- layer assumptions are 
applied in Cartesian coordinates, Pr t is given by 


Pr = ~~ c 
t V 


I 1 + f*t Cl ' C 62 1 ly l 


(2.81) 


where 


♦t = c ei + W + 

f = 0.41 f V 1/2 * x' 1 / 2 ) 2 

The last term in the brackets in the expression for f 1 is to include the 
effect of the step face in rapid expansion flow. For normal boundary- 
layer flow, this term would be deleted [137]. 

D. Engineering Parameters 

This section introduces the form of the engineering parameters used 


later in this study. 
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It is possible to derive expressions for Pr in a method similar to 
that used to derive an expression for c^. First, transport equations 
for the three turbulent scalar fluxes u^T* are derived. The gradient 
terms in these equations are removed by modeling assumptions to leave an 
algebraic expression for Pr^. libson and Launder [136] assumed that the 
transport of the scalar fluxes was proportional to the transport of k as 
Rodi did to derive an expression for c^. 

Watkins and Gooray [115] derived an expression for Pr in 
streamline and Cartesian coordinates. Their derivation was for a 
general elliptic problem. When the boundary- layer assumptions are 
applied in Cartesian coordinates, Pr t is given by 


n _ 1 ..A ,, .3u, 

Pr t = i;y i W 1 - c 92>^i 


(2.81) 


where 


♦t = C 81 + W + 


£' 


0.41 


^V 1 ' 2 + x’ 1 ' 2 ) 2 


The last term in the brackets in the expression for f' is to include the 
effect of the step face in rapid expansion flow. For normal boundary- 
layer flow, this term would be deleted [137]. 


D. Engineering Parameters 


This section introduces the form of the engineering parameters used 


later in this study. 
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1. Wall shear stress 

The shear stress at the wall for internal flow can be evaluated in 
two ways. The first way is by using the definition of t for a Newtonian 
fluid as follows 


T 

w 


(x) = T (x , 0) = u|^(x,0) 


(2.82) 


For internal flows, a second expression for can be derived by 
equating the forces and momentum fluxes entering and leaving a control 
volume spanning the channel width. The following expression results 


V x) 


. -m dp . m dM. 
(2 r -r* + r -T-) 
o dx o dx 


(2.83) 


where 

r 

M = / °prVdy (2.84) 

0 

Equation (2.83) applies since, for the boundary- layer equations, the 
pressure gradient is aligned with the x-coordinate and so can be 
expressed as a total derivative. The fact that the two expressions for 
T w should be equal, was used as an internal check in the computer 
program. 


2. Bulk temperature 

The mean or bulk temperature for variable property flow is given by 
the equation 

r r 

T, = (/ °pr m uc Tdy)/(/ °pr m uc dy) 

b 0 P o P 


(2.85) 


81 


A second expression for the bulk temperature can be derived by an energy 
balance on a control volume spanning the channel. The resulting 
expression is 


W = x' \ dx • 2 \ 0p 2 U 2 dy + I \ Vl dy + (2 86) 

r r 

V*i) 1 °P 1 u ( x 1 »y) c pl d y / J °P 1 u( x 1 »y) c pl d y 


The "l" and "2" subscripts refer to two different x-positions. The two 
expressions for T^ should be equal and so provide an internal check on 
the present predictions. 
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III. METHOD OF SOLUTION 

This chapter describes the method used to solve the equations 
developed in the previous chapter. The computational grid will be 
presented. The finite-difference discretization will then be discussed 
followed by the method used to solve the resulting system of equations. 
The method used to discretize the boundary conditions will also be 
presented. Finally, convergence, truncation error, and stability will 
be discussed. 

The equations are discretized iij.an implicit manner. When using 
implicit methods for parabolic marching problems, a system of equations 
must be solved for the unknowns at the next station beyond the known 
values. It is generally felt that implicit procedures are well-suited 
for parabolic marching problems [4] . 

A. Computational Grid 

The finite-difference equations used arc valid for uniform and 
nonuniform grids. A representative example of the orthogonal Cartesian 
grid used for this study is shown in Fig. 3. The mesh shown in Fig. 3 
is much coarser than those actually used. Grid refinement studies 
showed that 81 to 121 y-grid points and 135 x-grid points were adequate 
to resolve the flow field. The j index is used to specify y-position 
with j = 1 corresponding to the points on the wall, and j = NJ 
corresponding to the points on the centerline. The i index is used to 
specify x-position. For x = 0, i = 1. 



Figure 3. Finite-difference grid 
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Constant y-grid spacing was used for laminar flows. When 
predicting turbulent flows, it was necessary to use variable y-grid 
spacing in order to have a small Ay near the wall and a larger Ay away 
from the wall. For turbulent flow, inverses of the general stretching 
transformations of Roberts [138] as given in Anderson et al . [4] 
provided the y-grid spacing. The exact form of the : tretching 
transformations is given in Appendix F. The y-grid spacing was fine 
enough near the wall to ensure that at least two to three grid points 
were in the laminar sublayer. 

A geometric progression x-grid spacing was used that was defined by 


Ax + = KAx 


(3.1) 


where K is a variable greater than or equal to one. Ax_, and Ax + are 
shown in Fig. 3. The distance between x = 0 and the first solution 
station AXj is given by 


Ax 


1 


Ax 


1 


stop 
NSTP 1 


K = 


x „ (1-K) 

stop 

1 _ k nstp * 


1 

K > 1 


(3.2) 


where NSTP is the number of steps to be taken in the x-direction, and 

X ^ is the largest x-value in the solution domain. For flow with 
stop 

reattachment, K was set to one in order to provide adequate resolution 
at the point of reattachment. A value of K not equal to one was only 


used in flows with no step. 
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B. Velocity/Stream Function (u-^) Variable Equations 

The stream function and the st-reamwise velocity were used as the 
hydrodynamic variables in most of the computations done in this study. 
This section describes the finite-difference discretization of the 
governing PDEs with u, H, k and c as the variables. 

The continuity and momentum equations were solved in a coupled 
manner in the present study. Since they are coupled, solutions for both 
u and tfi (or u and v for the u-v equation set) are obtained 
simultaneously. This is different from the usual method of solving the 
boundary- layer equations. Usually, the momentum equation is solved 
first to find u using lagged values of v from the previous marching 
station. Then, the continuity equation is solved using the recently 
obtained u values to find v. By using the coupled procedure, the v (or 
tp) values are not lagged in the momentum equation so the predictions of 
u will be more accurate. Solving the momentum and continuity equations 
in a coupled implicit way is more complicated since the resulting system 
of equations is block tridiagonal with the tridiagonal elements bein' 
two by two matrices. When using an uncoupled implicit method, a 
tridiagonal system results where the tridiagonal elements are only 
single elements. Kwon and Pletcher [39] reported that the velocity 
profiles and pressure gradient predicted by the uncoupled scheme showed 
wiggles in regions of flow separation. When they coupled the continuity 
and momentum equations, the wiggles disappeared. 
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The energy, k and e equations were solved uncoupled. The 
continuity and momentum equations were solved first to obtain u and at 
a specific x. If the temperature f-eld was desired, or if the flow was 
turbulent, the u's and ^'s were used in the energy, k and z equations. 


1. Continuity equation finite -difference discretization 

The discretization of the continuity equation will be discussed 
first because it does not have the same general form as the other four 
PDEs. Figure 4a shows the difference molecule used for the continuity 
equation, Eq. (2.14). The finite-difference approximation of Eq. (2 j. 4) 
is 


, m .i+1 _ 1 , . i+1 .i+1. 

Cpr u Vi = j ' ♦j-i ) 


(3.3) 


When a subscript or a superscript is i or j plus or minus 1/2 as in Eq. 
(3.3), the value of the variable is taken as the average of the values 
at the two adjacent nodes, i.e., 

Vi - (0 j + w 2 - 


The finite-difference form of the continuity equation is then 


1, i+1 m i+1 i+1 m i+1. 3. , , i+- 1 .i+1. 

~(p. r.u. + p. _r. .u. ,) = - — (^. - V*. ,) 

2 J J J J’l j-i j-i &y_ 3 j-i 


(3.4) 


At this time, the density, p, is unknown at the i+1 station. In order 
to decouple the momentum and continuity equations from the energy 
equation, the density values will be lagged in the x-direction by 


setting 
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i+1 


1 

P J 


(3.5) 


The fluid properties at the i+1 marching station will always be lagged 
in this way. The continuity equation is now of the form 


, i+1 , i+1 i+1 i+1 . 

b.u. , + e.u. - 9 . + 0. , = 0 

J J-l J J J J-l 


(3.6) 


The expressions for b^ and e^ are found in Appendix B. 


2. General finite -difference discretization 

The momentum (Eq. 2.20), energy (Eq. 2.27), k (Eq. 2.68), and t 
(Eq. 2.68) equations can all by written in the same form as Eq. (2.68) 
as follows 


, ,1$ 1 30 30 - I 3_ f r 30, . q 

Cp U 3x r 3x 3y r 3y [ F 0 3y ] S < 


(3.7) 


The values of the T 's and the S.'s for each of the different equations 

0 0 

are given in Table 4. It is issunced at this point that y^ and Pr^ have 
been found through the appropriate turbulence model as discussed in 
Section II.C. The m superscripts on the r's have been dropped since it 
is more efficient to simply set r(y) = 1 for planar geometries and let 
r(y) be the distance from the centerline for axisymmetric geometries. 

The terms enclosed by the dotted "ertical lines in the k and e equations 
shown in Table 4 are the low-Reynolds-number terms added by Chien [135] 
to make the model valid in regions having a low Reynolds number of 


turbulence. 


89 


This section discusses those points of the finite-difference 
procedure that are common to all of the equations (momentum, energy, k, 
and e) fitting the form of Eq. (3.7). In the following sections, the 
details specific to a particular equation will be discussed. 


Table 4. 


Expressions for S. and IV 

♦ 9 


for the general equation 


* 

% 

s * 

u 

u t + 

-dE 

dx 

H 

H- + 

Pr Pr 


k 

*■ + “t 

f iu N 2 i 2uk i 

*tV 

e 


c £ u . _£ r c > f > E +I 2yk e (-c 4 y + )i, 

1 k W t C 2y J p k l 2 ,f • 'py* ,J 



« - 1 - 


Six different finite-difference molecules were used to discretize 
Eq. (3.7). The particular molecule used depended on the direction of 
the flow. Figure 4 shows the six different molecules and the direction 
of the "wind" for each. The molecules shown in Figs. 4b and 4e which 
use central differencing in the y-direction were used for the majority 
of the grid points. The others were used only when the velocity in the 
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y direction caused nonviscous behavior when the molecules of Figs. 4b or 
4e were used. This nonviscous behavior is discussed further in the 
Stability section. 

a. First global iteration using FLARE Equation (3.7) in 
finite-difference form for the case when the FLARE approximation is in 
effect or when the flow is in the positive x-direction is 
1 i i+1,, i+1 ,i, 1 , . i+1 . i^ r *,i+l _ 

4r cp j u j ( *s ' V ' rsr ( *j * *s >Vj " 

S + i- D(r ,r »‘ +1 ) (3 - 8) 

♦ r j J J J 

where 6^ is a central difference operator valid for nonuniform grids 
given by 


6 % = -a ) + -4 )| 

Vj Ay + +Ay_Uy + ^j+l V Ay_^j *j-l ; | 


(3.9) 


and D is a specialized diffusion difference operator used to discretize 
the term 


— ( r r — ) 


The form of D is 


D(r rVi +I) = 


(3.10) 


? / Ay_._ p , i+1 i+1 + 

A y+ +A y. \(Ay + ) 2 j + ^ 9 j + l V j Uy_ Ay + JAy + +Ay_ 

Uy + l, j+l ’j } Ay_ lf j (AyJ 2 r j-r j-r’j ? j-l ; / 
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b. After first global iteration After the first global 
iteration, the FLARE approximation was discarded and the x-convective 
term of Eq. (3.7) was differenced in a local upwind manner according to 
the local streamwise velocity. If u was positive then the difference 
molecule shown in Fig. 4b was used; if u was negative, the molecule 
shown in Fig. 4e was used. For the case with negative u, the finite- 
difference discretization of Eq. (3.7) is 


*»+ j j j j ' 


1 ✓ ,i+2 . i+1..* i+1 

r Ax y v j 


s . + 


(3.11) 


This allows downstream flow information in the global sense to influence 
the solution in the reversed flow region. For positive u, the 
discretization procedure was the same as that used for the first global 
iteration with FLARE. 

For turbulent flow calculations, the Newton linearization was 
dropped after the first global iteration and the nonlinear terms were 
linearized using values from the previous global iteration. 

c. Upwind discretization Equations (3.8) and (3.11) can both 
be written in the form 


, .i+1 , , .i+1 , A i+1 

Vj-i + Vi + “jVi 


= c , 


(3.12) 


For positive d ^ , the coefficients b^ and a^ must be negative for Eq. 
(3.12) to correctly model viscous behavior [4]. For laminar flows, 
nonviscous behavior was not a problem. However, for turbulent flows, a. 


• 4 - 


92 


and bj did become positive at positions near the step and near the 
reattachment point. 

Upstream differencing of the y-derivative in the y-convective term 

when either a. or b . was positive ensured that both a. and b. were both 
J J J J 

negative. Central differencing was still used in the diffusion operator 
D. If a^ was positive, the finite-difference molecules u c H and "f" of 
Fig. 4 were used; if b. was positive, molecules "d M and "g" were used 
for the y-convective term. This is easily done by substituting the 

# 

correct upwind finite-difference for the central difference operator 6^ 
as 


,* . i+1 


. i+1 „ i+1. 


6 => (0. -0. , )/Ay , for a. > 0 

y J J j-i - J 

* i+1 -v ^ i+1 ^ i+1 N/A 

S y*j => )/4y + ’ 


for b . > 0 
J 


3. Momentum equation finite-difference discretization 

In this study, the momentum and continuity equations were solved to 
the wall for laminar and turbulent flows. Most other investigators used 
wall functions near the wall for turbulent flows so that the solution 
point nearest the wall is well-away from wall effects. A turbulence 
model valid in reversed flow that is accurate near the wall as well as 
in the high-Reynolds -number region away from the wal highly 
desirable. 

a. Linearization For the momentum equation with 0 = u, Eq. 
(3.8) is nonlinear since two unknown u values at the i+1 station are 
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multiplied together. There are several ways to linearize the 
coefficients of Eq. (3.8): 

• lag the coefficients 

• extrapolate the coefficients 

• update of the coefficients by simple iteration 

• update of the coefficients by Newton linearization 

A detailed description of each type of coefficient linearization can be 
found in Anderson et al. [4]. 

Kwon and Pletcher [39] studied the effect of several linearization 
schemes. They reported that only Newton linearization while coupling 
the continuity and momentum equations gave well-behaved predictions when 
there were large areas of recirculation. Newton linearization also 
greatly enhanced the rate of convergence of the linearization iteration. 
Due to the findings of Kwon and Pletcher, Newton linearization with 
coupling was the linearization scheme used in the present study. 

The main idea behind Newton linearization is to replace the 
coefficients of the convective terms that cause the nonlinearity with 


i+l = ji+l + 5 
J J ♦ 


(3.13) 


where 6 is a small change between the converged value at the i+1 level 

and is a provisional value from a previous Newton linearization 

iteration. After each of the variable coefficients of Eq. (3.8) is 

2 

replaced with Eq. (3.13) and the 6^ terms dropped, a linear equation 
results. For the first Newton iteration, the provisional values are 
lagged and the solution is stepped from the i marching station to the 
i+1 marching station. The provisional values are then updated with the 
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predicted values at. the i+1 level and the equation is solved again for 

the unknowns at the i+1 marching station. In uhi's study, the iteration 

-4 

was continued until the average solution change was less than 5.0x10 
The average solution change was defined as 


maximum{Z. |u* + *-u. |/Z. 

J J J J 


j J J J J 


The solution was then advanced to the i+2 marching station in the same 
way. 

After using Newton linearization to linearize Eq. (3.8) and 
grouping the coefficients for each of the different unknowns, the 
following equation results 


_ i+1 . n i+1 , . i+1 . r .i+1 

B j u j-i + D j u j + Vj« + E /i 


H.X + C. 
J J 


(3.14) 


The pressure gradient, -dp/dx, has been expressed as x. The 
coefficients of Eq. (3.14) are given in Appendix B. 

After the first global iteration, if u was negative, Eq. (3.11) was 
used ip place of Eq. (3.8). Equation (3.11) was linearized in the same 
way as Eq. (3.8). After grouping the coefficients of the unknown 
variables, the same form as given by Eq. (3.14) results. The variables 
with i+2 superscripts a::e considered as knowns since the values used are 
obtained from the previous global iteration. The coefficients of Eq. 
(3.14) for the case of negative u are given in Appendix B. 

b. Boundary conditions The boundary conditions for u and 4 > for 
use in solving the momentum and continuity equations follow. At the 
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wall (y = 0), Uj ^ and ^ ^ are simply set to zero. At the 
centerline, 3u/3y = 0 is specified in the following way. The equation 
valid on the centerline for both axisymmetric and planar geometries is 
(see Eq. 2.20 and 2.21) 


3u „ . 3 . . ,3u 

CpU ^ “ * + 2 3y (U+ Va? 


(3.15) 


To specify the symmetry condition at the wall, a reflection boundary 
condition is used by defining a pseudopoint, y^j + ^> beyond the y-value 
corresponding to the centerline, y^j- Equation (3.15) is finite 
differenced using the molecule of Fig. 4b to give 


1 i i+1, i+1 l v . ~m„,. _ i+1 

Ax_ p NJ U NJ U NJ " U NJ ) " X 2 D(1,r u,NJ ,U NJ ) 


i+1 


(3.16) 


Due to symmetry, for j = NJ 


u 


( ) _ .,( ) 


NJ+1 


" U NJ-1’ Ay - “ Ay +' r u,NJ+i ~ r u,NJ-i 


(3.17) 


Substituting Eq. (3.17) into Eq. (3.16) gives 


®NJ U NJ-1 + D NJ d NJ “ + C NJ 


(3.18) 


The expressions for the constants cf Eq. (3.18) are given in Appendix B. 


4. Solution of coupled f inite-dif ference hydrodynamic equation system 
This section outlines the method used to solve the finite- 
difference approximations of the continuity and momentum boundary- layer 
equations for compressible variable property flow with u-^> as the 
variables. Since Eqs. (3.6) and (3.14) both involve uj*|, u* +1 , 
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and ^ , they must be solved simultaneously. For j = 2 to NJ-1, the 

following finite-difference equation holds 


B. 0 
J 


b. 

LJ 


i+1 

U j-1 


„ i+1 

\j-i 


> + 


D. 

J 


. e . -1 

L J 


i+1 

u. 

J 


i+1 


'j 


> + 


A. 

J 


o ' 1 


0 


i+1 

\»+l 


.i+1 

♦j+1 


> = < 


H.X+C. 
J J 


(3.19) 


Assembling the set of equations consisting of Eq. (3.19) written for 
each y-grid point at a given x-position results in a system of linear 
equations that must be solved simultaneously. The resulting system is 
block tridiagonal with each block consisting of a two by two matrix as 
follows 


(3.20) 



where 
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• * 
A, 0 


» * 
B. 0 


D j E j' 

II 

j 

; [BJ j * 

j 

; [BJ j - 


1 c 

o 

b j *. 


e . -1 

J J 


(Oj = 


[H.X + C. 
J J 

; (u>. = 

l+li 

“j 

G.X + d. 

j r 

j 

.1+1 
*J J 


(3.21) 


The coefficients of Eq. (3.21) are the same as those listed in Appendix 
B. 

Equation (3.20) was reduced to upper triangular form by using a 
modified Thomas algorithm. The Thomas algorithm is commonly used to 
efficiently solve tridiagonal systems of equations in which the elements 
are single coefficients. The same method can easily be applied to solve 
Eq. (3.20) due to the sparseness of the submatrices [A]^ and [BJ^. To 
eliminate the submatrices below the diagonal on the jth row of 
submatrices, the j-1 row is multiplied by *[B]j[D]^^ and added to 
the jth row. [D]^ the * nverse °f [ D ] ^ . Since it is known that 
[B]^ will be zeroed and that [A]^ will not change because all the 
elements above it are zero, it is only necessary to modify [D]^ and 
(C}j. The elimination process is started with j = 2 and proceeds to j = 
NJ. The new diagonal and right hand submatrices become 


IDlj = + [DJj 

-1 


(3.22) 

(C)" = -(BljJDljl^Oj.i + {C}j (3.23) 

The exact form of the new submatrices is given in Appendix C. 
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The pressure gradient parameter X has been grouped with known 
coefficients on the right-hand side of Eq. (3.20) even though it has 
been unknown to this point. However, after the system of equations is 
in upper triangular form, the mass flew constraint which specifies at 
the centerline (Eq. 2.30) is implemented to solve for the unknown 
pressure gradient, X. This is done by noting that after Eq. (3.20) is 
reduced to upper triangular form, the reduced form of the momentum and 
continuity equations for j = NJ is 


‘“InjKj 1 ' < C >NJ 


(3.24) 


★ if 1 

(C) N j contains the unknown X; {U}^j contains only the unknown u^j. 
since is a known boundary condition. Since there are only two 
unknowns and two equations, u^** and X can be algebraically 
determined. The resulting expressions for u^j^ and x are given in 
Appendix C. 

After determining u^ 1 and X, a matrix form of back substitution 
was used to find {U}^ for j t NJ as follows 


(U >j - U»lj)' 1 «C)‘ +1 - (A)j(U) J+1 ) 


(3.25) 


Appendix C shows the specific form of the back substitution. 


5. Energy equation finite -difference discretization 

As was done when discretizing the momentum equation, two types of 
discretizations were used depending on whether the FLARE approximation 
was being used and the sign of u. After the first global iteration, the 
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FLARE approximation for the energy equation was discarded. The x- 
derivatives in the energy equation were then differenced in a local 
upwind manner to allow energy to be convected from the globally 
downstream direction when flow reversal was present. Upstream 
differencing for the y-derivative in the y-convective term was used as 
necessary to ensure viscous behavior of the predictions. 

The discretisation of the energy equation for a general point 
removed from the boundaries is given by Eqs . (3.8) and (3.11) with 
♦ = H. Equation (3.8) was used for the first global sweep and when u 
was positive; Eq. (3.11) was used when for the second and subsequent 
global sweens when u was negative. The u's and if/'s are known since the 
momentum and continuity equations were solved first. Thus, the 
resulting finite-difference equation is linear in H. The source term of 
the energy equation as given in Table 4 is discretized as 



= — D(r.,S.,-Hu i+ V) 

J J J 


When the coefficients of the unknown Hs are grouped, che result is 


, „i+l , , „i+l , „i+l 
b.H, , + d.H. + a.H... = c. 
j J-l J J J J+l J 


(3.26) 


The expressions for the coefficients are given in Appendix D. 

a. Boundary conditions Now that the governing FDEs away from 
the boundary have been determined; the appropriate boundary conditions 
must be discretized in order to obtain the complete system of governing 
equations. The requirement that the flow be symmetric about the 


100 


centerline is used as the boundary condition at the centerline. The two 
possible boundary conditions at y = 0 are: a specified surface 
temperature (Eq. 2. 32a), and a specified heat flux (Eq. 2.32b). 

At the centerline, the symmetry condition leads to the following 
equation for j = NJ 


1 i i+l,„i+l 
Ax_ p NJ U NJ (H NJ 





(3.27) 


The symmetry condition states that 


Using this fact 


with Eq. (3.17), Eq. (3.27) becomes 


. „i+l . j H i+1 _ r 
b NJ H NJ-l d NJ H NJ NJ 


(3.28) 


The values of the coefficients are listed in Appendix D. 

1) Specified wall surface temperature For a specified 
temperature boundary condition at the wall, ^ = c^T^x) . The wall 
boundary condition is implemented by using Eq. (3.26) for j = 2 as 
follows 


, „i+l . „i+l 
d 2 H 2 a 2 H 3 


= c. 


- S H 1+1 
b 2 H l 


(3.29) 


Since is known, it has been moved to the right hand side with ca- 

using Eq. (3.29), Eq. (3.26) for j-3 to NJ-1, and Eq. (3.28) for j = NJ 
gives a tridiagonal system that can easily be solved by Gaussian 


elimination. 
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2) Specified wall heat flux For a specified heat flux at 
the wall, a one-sided three point finite-difference approximation of 
3T/3y in Eq. (2.32b) gives 


A U i+1 4. U i+1 

d l H l + “l H 2 


+XH 


i+1 


= c. 


(3.30) 


where d^, a^, X, and c^ are given in Appendix D. In obtaining Eq. 
(3.30), 0 ^ was assumed to be locally constant. 

Equations (3.30), (3.26) for j = 2 to NJ-1, and (3.28) for j = NJ 
make up the system to be solved for the case of a specified wall heat 
flux. Eq. (3.30) destroys the tridiagonal form of the system of 
equations. This is easily overcome by adding - X/a ^ times Eq. (3.26) 
with j = 2 to Eq. (3.30). This modification of the first row of the 
coefficient matrix and right hand column matrix does not harm the 
diagonal dominance of the algebraic system. The Thomas algorithm can 
now be used to efficiently solve for the unknown Hs at the i+1 x- 
position. 

b. Initial conditions At x = 0, initial conditions must be 
specified for H. For an unheated channel upstream of the step, the 
initial condition for the temperature above the step is simply 


T(0,y) = T 


inlet 


Since the experimental data sets to be compared with had an insulated 
step, the initial condition below the face is 3T/3x = 0.0. For the 
first sweep down the channel, 3T/3x in the x-convective term of the 
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energy equations was set to zero. This removed any influence of the 
temperature along the face of the step for the first sweep down the 
channel. For subsequent sweeps, a three-point finite-difference 
representation of 3T/3x =0.0 was used to extrapolate backwards to the 
face of the step from the first two x-grid stations beyond the step. 

The extrapolated face temperatures were then used in the next sweep down 
the channel. 


6. k-r equation f inite-dif ference discretization 

The e - equation was solved before the k-equation since it is the 
more approximate of the two equations [3] . Lagged values of k or values 
of k from a previous iteration were used in the source term of the e 
equation to uncouple the two equations. The k-equation was then solved 
using the recently computed values of £. Upon knowing k and c , y^ was 
then calculated by Eq. (2.67). 

The convective and diffusive terms of the k and e equations were 
finite differenced using Eqs. (3.8) and (3.11) with upwind differencing 
of the y-derivative in the y-convective term as needed. The resulting 
equation is of the form 


. , i+1 . . . i+1 , . i+1 

b.k + d.k. + a.k. , , 
3."1 J J J J+l 


= c , 


(3.31) 


The values of the coefficients in Eq. (3.31) and the similar equation 
for t are listed in Appendix D. 

a. Source terms Unlike the momentum and energy equations, the 
source terms of the k and e equations need special treatment. As the 
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wall is approached, the source terms begin to dominate the diffusion and 
convection terms [122]. Special handling of the source terms in the k- 
equation is needed to ensure that k remains positive [139]. 

The source term of a PDE must be written in the following linear 
form when approximating the PDE: 


S. = h * + S * , 
♦ $,d 9,c 


(3.32) 


where S. , is the part that is included in the d coefficient of Eq. 

9 jd 

(3.31) and S is the part that is included in c . of Eq. (3.31). If 
v > o j 

S. is always positive and S. , is always negative, for initial 
♦ ,c $,a 

positive 0 will always be positive. The source term of the k- 
equation written in the form of Eq. (3.32) is 


S R = (“H)k + (2p t (|^) 2 - pt) (3.33) 

If ^ was positive, Eq. (3.33) was used as the expression for the 
source term. When S, became negative, S was set to zero and S , 

K. y C K | C K j O 

took the form 


k,d 


■p + (2 vfp>’ - 


(3.34) 


ir 

where k is a lagged value of k or a k given by a previous iteration. 
The c -equation source terra is 


S 

E 



* 

fe 


2 y,. , / 3u . j 

p )E + VtV 


* 

E 

* 


k 


(3.35) 
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Equation (3.35) shows that S £ c is always positive and S £ ^ is always 
negative as they should be to ensure positive c. 

The previous discussion of the source terms for the k and t 
equations included the terms added by Chien (135) for low-Reynolds- 
number regions. 7f the low-Reynolds -number form is not desired, the 
terms in the vertical dotted lines of Table 3 must be deleted. 

b. Near -wall models For the high-Reynolds-number k-e model 
that was used for all the separated flow calculations, the near-wall 
region must be modeled differently than the rest of the flow. There 
were three different methods used for specifying the turbulent viscosity 
near the wall: (1) the maximum-shear-stress model, (2) the Prandtl- 
mixing-length model and, (3) the Prandtl -mixing- length model with 
variable A + . 

The maximum shear stress model is based on the inner viscosity 

model of Johnson and King (121). The point at which the switch was made 

from the wall model to the k-c high-Reynolds-number PDE model was at a 

constant y value, y^. The value of y^ was specified as the y 

corresponding to a y + of 30.0 to 50.0 for fully developed turbulent flow 

at the same Reynolds number (based on D and and the average outlet 

velocity). The predictions were not sensitive to the value of y^ within 

this range. The turbulent viscosity in the near wall region for method 

1 was given by Eq. (2.65). Equations (2.76) and (2.77) give the 

expressions used for k(y, ) and c(y, ). 

d b 
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For method 2, the boundary conditions were applied at a y + value of 
30.0 to 50.0 with similar results for both values. For method 3, the 
boundary conditions for k and e were applied at 

y + = 2.0 A + 

For method 3, A + was varied according to Eq. (2.64). In the 
recirculation region and near reattachment, A + reached values as low as 
1.0. Ordinarily, this is well-into the viscous sublayer and much too 
small a value of y + for the high-Reynolds-number k-e equations to be 
valid. Johnson and Launder [113] found that by reducing the thickness 
of the laminar sublayer below that normally found in turbulent boundary 
layers, the wall functions of Chieng and Launder [111] gave much better 
predictions of the heat transfer in recirculating flow. 

Since the value of y at which the k-E boundary conditions were 
applied varied for methods 2 and 3, it was necessary to specify values 
for k and r. for y values corresponding to y + less than 30. For methods 
2 and 3, the k and t values for y + < 30.0 were given by Eqs. (2.74) and 
(2.75). Equation (2.69) then gave the near-wall turbulent viscosity. 

c. Initial conditions To start the marching procedure for the 
first sweep down the channel, lagged values for k must be used in the 
source terms. Since the source terms of the k and e equations contain 
k \ zero k values are inappropriate. To remove the singularity, the 
values for k and z on the face of the step were set to the values of the 
point just above the lip of the step at the inlet to the channel (j = 
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NH+1, i = 1 as in Fig. 3). Sometimes these values caused the solution 
to diverge for fine grids and large expansions as the solution was 
marched from the step on the first iteration. This divergence could 
sometimes be overcome by using the k and t values at a point further 
above the lip for the face values ( j = NH+2 or NH+3) . 

C. Primitive Variable Hydrodynamic Equations 

This section describes the method used to solve the constant 
property continuity and momentum equations in a coupled manner using 
primitive variables, that is, without introducing the stream function. 
Previous primitive variable boundary- layer calculations for separated 
flow with the momentum and continuity equations uncoupled predicted a 
solution with small unphysical "wiggles" {39]. The main reason for 
predicting the flow with the primitive variable equations was to see if 
the use of tp was essential to obtain a satisfactory solution when large 
regions of separation were present. It will be shown in Section 
IV.A.l.d that the primitive variable formulation gives predictions 
identical to those of the u-^ predictions under similar assumptions. 
This indicates that the equation coupling overcomes the small 
oscillations and not the choice of variables. Only the planar two- 
dimensional equations for constant property laminar flow will be 
developed. The predictions of the primitive variable formulation would 
give the same predictions as the u-i> variable case for axisymmetric 
geometries, for problems with heat transfer, and for variable property 
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flows. The effect of iterating globally over the flow field was not 
studied for the primitive variable case. 

For planar 2-D coordinates, the transverse dimensionless velocity V 
can be given in terms of the dimensionless stream function as 


V = — = v ^ e 

~3X u. 

1 


Substituting Eq. (3.36) into Eq. (2.45) gives 


rII 3U 3U = i_ r au. 

U 3X V 3Y P 3Y ( 3Y ) 


(3.36) 


(3.37) 


The continuity equation is 


3U 3V 
3X 3Y 


0 


(3.38) 


Equations (3.37) and (3.38) are for constant property laminar flow and 
are nondimens ionalized in a way so that they are independent of the 
Reynolds number. 


1. Finite -difference discretization 

Equation (3.38) is finite differenced in the following manner 


Hi i (u j +i - u i> + (u £ - u }-i» + v r - v ?-i = ° 


,i+l 


i+1 


,i+l 


(3.39) 


The finite-difference discretization of Eq. (3.37) is 


1 -rri+l/rri+l „iv , 1 ,,i+l ,„i+l ,,i+l. 

4X _0U. (U. - v ♦ (U. +1 - Uj.j) = 


p + 


1 -(uft! - uh 


AY + +AY_ AY + j+1 j' 


AY 


1 (U 1+1 - U 1 + J)] 

J J-l 1 


(3.40) 
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Note that the FLARE approximation is in effect to allow the solution to 
be marched through regions of reversed flow. Equation (3.40) is 
linearized with Newton linearization as was done for the u-1> variable 
scheme. Equation (3.39) and the linearized form of Eq. (3.40) can be 
combined into a coupled equation as 


(BI j {u) j-i + (D] j (u> j + [A V u) j+i = {c} j 

where 


• ■ 
A. 0 


B. o' 


"d. 

8.1 

J 

; (BJj = 

J 

* [D1 j - 

j 

J 

0 0 
• « 

b. -1 
L J J 

L e j 1 J 



H.B + C. 
J J 

G.B + d. 
• J J 



i+1 
u . 



(3.41) 


(3-42) 


Writing Eq. (3.41) for each value of j gives a linear system of 
equations. Appendix E gives the fen? of the coefficients of Eq. (3.42). 


2. Solution of the system 

Two different methods were used to solve the system of equations 
consisting of Eq. (3.41) for different values of j. Regardless of the 
method used, the pressure derivative, 8, was obtained as part of the 
solution similar to what is done when using an inverse boundary- layer 
method. Two different methods were used to solve for 8: solving for 8 
algebraically as was done for the u-^ variable method, and solving for 8 
using an iterative secant numerical procedure. 
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a. Modified Thomas algorithm When Eq. (3.41) is assembled for 
all values of j, the resulting equation is the same form as Eq. (3.20), 
and so can be solved by the modified Thomas algorithm. The subraatrices 
below the diagonal are eliminated by row operations on each row of 
submatrices. Appendix E gives the resulting values of the components of 
the submatrices. 

After the system of equations is in upper triangular form, the 
coupled momentum and continuity equations for j = NJ are given by 



SnjIKj 1 ] 


["NJ 6 + C Nj| 


g nj 6 + d nj! 


The primes denote that the coefficients are for the upper triangular 

system of equations. Since V*.** = 0, Eq. (3.43) has two unknowns, 3 

NJ 

and which can be obtained algebraically. Now a matrix form of 

back substitution (Eq. 3.25) is used to solve for the unknown Us and Vs 
at the i+1 marching station. Appendix E gives the exact form of the 
submatrices used in the back substitution. 

b. Pressure -derivative secant If 3 were known, the system of 
equations consisting of Eq. (3.40) for each grid node could be directly 
solved by any suitable solution scheme for a system of equations. To 
start the secant procedure, 3 is simply guessed and the system of 
equations solved. If the guessed 3 is correct, then the boundary 
condition that = 0 will be satisfied. For each incorrect 3 guessed, 
there will correspond a nonzero V^j. After two guesses for 3 , a secaiit 
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procedure can be used to predict a new guess for $ that will yield a 
value of V VIT closer to zero. 

Three different variations of the secant method were tried: 

1. A secant procedure to find the correct pressure derivative 
was nested in the Newton linearization loop. 

2. The Newton linearization loop was nested in the secant loop 
to find the correct pressure derivative. 

3. The Newton linearization loor was removed resulting in lagged 
coefficient linearization with a pressure secant procedure to 
find the correct pressure derivative. 

Figure 5 shows a block diagram of the method used for variation 1. 
For variation 1 above, the provisional values used in the Newton 
linearization were set to the values from :he previous channel station. 

A secant method was used to determine the correct value of the pressure 
derivative in the following way. Two guesses for the pressure 
derivative were obtained by multiplying the pressure derivative at the 
previous station by 1.04 and 0.96. The pressure gradient for the first 
step was taken as the pressure gradient that would occur in fully 
developed channel flow. The general block solver NBTRIP provided by 
Sukumar R. Chakravarthy of Stanford University [4] provided two 
solutir,.s to the system resulting from these guessed pressure gradients. 
For the correct pressure derivative, the transverse velocity at the 
centerline, V., T should be zero. The centerline Vs obtained from the 
first two guesses were used in a secant procedure to predict a new 
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Figure 5. Block diagram for the pressure secant algorithms 
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pressure derivative. NBTRIP gave a new solution based upon this new 
pressure derivative. The secant iteration continued trying to force V 

IN J 

to zero until the relative change of the pressure derivative was less 
than some small prescribed value. At this point, the provisional 
variables used in Newton linearization were updated and the process 
repeated. Thus, tie pressure iteration was nested in the Newton 
linearization iteration for each step down the channel. 

Variation 2 is similar to variation one except that the Newton 
linearization loop was nested in the pressure secant loop. A first 
guess for the pressure derivative was chosen in the same way as in 
variation 1. NBTRIP gave a prediction for the velocity profiles at the 
new step from the input guessed pressure derivative and the velocity 
profiles at the previous step. The Newton linearization algorithm then 
looped until the provisional velocities stopped changing. At this 
point, a second guess for the pressure derivative was input to the 
linearization loop and the process repeated. The secant method then 
predicted a new pressure derivative. This new pressure derivative was 
used as input to NBTRIP to predict new velocity profiles and the 
linearization iteration was repeated. The pressure loop (with 
linearization loop nested in it) was repeated until the pressure 
derivative change was very small. 

Variation 3 only had the pressure secant loop. It is the same as 
variation 1 except that there was no Newton linearization loop (see Fig. 


5 ). 
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Section IV.A.l.d gives the results of the different primitive 
variable schemes. 


D. Convergence 

With any numerical approximation to a partial differential equation 
(PDE), one must be concerned with the errors in the approximation. It 
is hoped that a finite-difference equation (FDE) solution is close to 
the exact solution of the PDE. If the FDE formulation is consistent and 
stable, it converges to the solution of the PDE as the grid is refined 

m. 

1. Consistency and truncation error 

A FDE formulation is consistent if the truncation error (TE’' goes 
to zero as the computational grid is refined. The truncation error is 
defined as 

TE = PDE - FDE (3.44) 

The truncation error was found by using Taylor series expansions of 

the derivatives in the PDEs. The truncation error of the continuity 
2 

equation is 0(Ay_); the truncation error of all the other FDEs is 
0(Ax,Ay Ay ,Ay -Ay ) when central differencing of 80/3y is used in the 

* T " r 

y-convective term. (The truncation error degrades to 0(Ax,Ay) for the 
upwind differencing of the y-convective term but this usually only takes 
place near the step where tne solution procedure is admittedly 
approximate.) Note that for uniform grid spacing, the TE becomes 
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2 

formally C(Ax,Ay ). Since Ax, Ay 4 0 as the mesh is refined, the FDE 
discretization is consistent. 

Even though the TE is only formally first order in Ay, it can be 
shown that the TE in the y-direc+ iovi behaves in a second order fashion 
as was suggested by Blottner [140]. The dominant term of the truncation 
error is 

«y + - iy.) ■ (y J+1 - * 1 i + yj.^fp <3 - 45) 

If y is defined by a stretching transformation, y = y(C), of the C space 
which has a uniform grid spacing, AC, then 

0 = AC* (y j+l ” 2y j + y j-l } + (3,46) 

Equation (3.46) is introduced to merely establish a link between the 
nonuniform y-grid spacing and the uniform C”£rid spacing. Equation 
(3.46) can be rearranged to give 

y j+ l ' 2y. + yj-1 = 0(AC 2 ) (3.47) 

The finite-difference discretization has a formal truncation error of 
2 

0(AC ) for the uniform grid space. Equation (3.47) indicates, that the 

largest term of the truncation error of the nonuniform grid behaves in a 

second order manner with respect to the uniformly divided C space. If 

more y-grid points are added, an additional corresponding number of 

C"grid points will be added. Therefore, the finite-difference 

discretization in the nonuniform grid space, y, behaves in a second 

2 

order fashion (TE = 0(AC )) as long as y can be defined as y = y(C). 
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2. Stability 

The second requirement for convergence is stability. A solution 
methc*d is stable if the round-off errors in the computer and the 
truncation errors do not grow as the solution progresses. For linear 
PDEs, a von Neumann stability analysis can determine if a FDE 
discretization is stabTe [4], The implicit method used in this study is 
unconditionally stable for linear initial-value PDEs and was assumed to 
be stable for the nonlinear boundary- layer equations of the present 
study. 

For linear initial-value PDEs that are solved in a marching manner. 
Lax's Equivalence Theorem shows that for a consistent FDE 
discretization, '.tability is a sufficient condition for convergence [4]. 
Lax's theo- em has not been proved for the nonlinear case but is 
generally accepted as valid. Since the present finite-difference 
formulation is consistent and stable, it is convergent. 

Even though most implicit marching schemes for parabolic PDEs can 
be shown to be unconditionally stable in the von Neumann sense, an 
inappropriate grid can cause wiggles in the solution due an unphysical 
modeling of the PDE. This is easily shown by looking at the following 
momentum equation for an uncoupled implicit marching scheme 

B.u 1+1 + D . u 1+ 1 + E.u 1+ J = C (3.48) 

J J-l J J J J+l j 

To realistically model viscous flows, B^ and should be opposite in 
sign to D. [4]. This gives the flowing mesh Reynolds number constraint 



v 


(3.49) 
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If Ay or Vj is too large, the mesh Reynolds number constraint is 
violated and unphysical behavior is predicted. 

Nonviscous behavior due to the violation of the mesh Reynolds 
number has not been documented for the case when the momentum and 
continuity equations are coupled. Kwon and Pletcher [39] have suggested 
that proper viscous modeling is achieved because v in the v3u/3y term is 
treated as a variable when the momentum equation is solved. 

The continuity equation is a mass balance and does not involve the 
viscosity so its effect when coupled with the momentum equation should 
not change the required character of the coefficients in the momentum 
equation. This suggests that and of Eqs. (3.19) and (3.31) 
should remain negative. 

Two common ways to ensure negative off-diagonal terms is to upwind 
difference the y-derivative in the y-convective term and to refine the 
grid. Since a large number of grid points were already being used in 
this study, upwind differencing was used to ensure proper viscous 
modeling. 
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IV. RESULTS 

This chapter presents the results of this study. The laminar 
results will be discussed first followed by the turbulent results. The 
laminar and turbulent sections are each divided into hydrodynamic and 
heat transfer sections. 


A. Laminar 

The laminar results were promising in general and shed light on the 
nature of the boundary- layer equations as compared to the Navier- Stokes 
equations. The solution algorithm is quite robust and well-behaved for 
all the cases studied. No artificial damping or overrelaxation of any 
kind was needed for the laminar predictions. A typical laminar run 
required 3-5 seconds of NAS AS/6 computer time per global sweep through 
the flow field. 

1. Hydrodynamic constant temperature 

The constant temperature laminar results are presented in this 
section. The continuity and momentum equations (Eqs. 2.44 and 2.45) 
were nondimens ionalized so that the Reynolds number dependence was 
removed. Calculations were made for expansion ratios of 1:1.2 through 
1:4 (d:D as in Fig. 2) for both symmetrical planar and axisymmetrical 
expansion flows. This range of expansion ratios included the range of 
those studied by others with experimental tests and numerical 
predictions for this type of geometry. Comparisons were made with 
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available experimental results and ..olutions of the Navier-Stokes 

equations. Particular attention was given to the influence of expansion 

ratio and Reynolds number on the level of agreement between the 

predictions of the boundary- layer equations and measurements or Navier- 

Stokes solutions. Some flow parameters were predicted better than 

others. Those generally we 11 -predicted were velocity and reattachment 

length. The higher the Reynolds number, the better the predictions. 

a. Convergence When marching through a region of reversed flow 

using the boundary-layer equations, it is important to ensure that the 

solution does not become unstable as the grid is refined. Blottner 

[141] reported that his solution to the s iender- channel equations 

diverged as he refined the grid when in the reversed flow region. For 

this reason, grid refinement was carefully studied to verify that the 

solution by the present method did not become unstable when the grid was 

refined. Figure 6 shows the effect of reducing AX on the reattachment 

length and the location and value of t . . Figure 6 shows the solution 

mm 

asymptotically approaches definite values as the grid is refined and 
does not diverge. 

Figure 7 shows the relative change of key parameters per global 

iteration. In Fig. 7, Z is a surrogate variable that takes on the value 

of X , the X-distance to f , , the Y-distance to T . , or T . . Since 
r min mm mm 

the relative change of these four parameters goes to zero with 
increasing iteration number (i), it is obvious that the global iteration 
converges to steady-state values. 




Figure 6. The effect of grid refinement on reattachment length, minimum stream function (f 
and distance to t for an axisymmetric 1:2 expansion 
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a-!)Z/ca-nz-a)Z] 


Iteration (i) 

Figure 7. Relative change of hydrodynamic flc parameters with global iteration for a planar 
channel expansion (Z is a surrogate iable) 
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b. Comparison with other data sets Reattachment length was 
we 11 ^predicted over a wide range of expansion ratios and Reynolds 
numbers. Figures 8 through 10 compare the reattachment lengths 
predicted by the boundary- layer equations with those of experiments and 
other mathematical models. The boundary- layer predictions of Figs. 8, 

9, and 10 are all for only one global sweep through the flow domain. 
Global iteration over the flow field had almost no affect on the 
reattachment length predicted by the boundary- layer method. 

Figure 8 compare-, the reattachment length and distance to the 
vortex center predicted by the boundary- layer equations using the FLARE 
approximation with the experimental and computational results of Macagno 
and Hung [5] for a 1:2 axisymmetric expansion. The comparison between 
reattachment lengths predicted by the boundary- layer solution and the 
results of Macagno and ilung are excellent for Reynolds numbers above 
twenty . 

Figure 9 compares the reattachment lengths from the boundary- layer 
solution with that of other investigators for different Reynolds numbers 
and expansion ratios for symmetric planar expansions. The reattachment 
length measured by Durst et al. [8] was taken from a photo of smoke in 
air from which it was difficult to tell accurately where reattachment 
occurred. From Fig. 9 it appears th-t the agreement is very good for 
h/d from 0.0 to 0.5 and at least marginal fror 0.5 to 1.0. 

Table 5 compares the predictions of the boundary- layer equations 
and the Navier-Stokes equations [43] for the location of flow 
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Figure 8. Comparison of the distances to flow reattachment and vortex center for experiment [5], 
Navier-Stokes solutions |5], and the boundary- layer equation solution with FLARE for a 
1:2 pipe expansion 
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Figure 9. Nondimens ional distance to reattachmer.t for a planar channel expansion as a function of 
expansion ratio predicted by the boundary- layer equations [37,27], the Navier-Stokes 
equations [26,43], and experiment [8] 
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reattachment, the location of the minimum stream function (Y . ) , and 
’ mm 

the absolute value of f , . The ratio of the inlet to outlet plate 

min 

spacing was 1:3; the Reynolds number was 37.3. Even for this large 
expansion ratio and low Reynolds number, the distance to reattachment 
predicted by the boundary- layer equations is within 5 % of that predicted 
by the Navmr-Stokes equations. 

Table 5. Comparison of Navitv-Stok.es and bour.dar- -’ayer predictions for 
a 1:3 planar expansion 


STUDY 

X 

r 

X tj T . ,Y to T . 

mm mm 

f . 
min 


Osswald et al. [43] 

0. 1030 

■j 0290,0.615 

0.0515 


. resent 

0.0981 

0.0227,0.647 

0.0668 



Figure 10 shows the axisymmetric analog of Fig. 9. Comparisons are 
made with Pollard's [12] and Macagno and Hung's [5] solutions of me 
Navier-Stokes equations. Again the agreement appears excellent for h/d 
less than 0.5 and good for values as high as 1.0. 

A fourth order polynomial fits the boundary- layer predictions of 
Figs. 9 and 10 to within 3%. The reattachment length can be expressed 
as 

X r = A + B(jj) + C(^) 2 + D(^) 3 + E(|j ) k (4.1) 

Table 6 gives the values of the constants of Eq. (4.1). Equation (4.1) 
is valid for 0.1 < h/d < 1.5. 
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Figure 10. Nondimens ional distance to reattachment for an axisymmetric expansion as a function of 

expansion ratio predicted by the boundary- layer equations and the Navier-Stokes 
equations [5,12,43] 
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Table 6. Coefficients for Eq. (4.1) 


Geometry A B C D E 


Symmetric Planar -0.001277 0.004262 0.1748 -0.09451 0.02086 

Axisyrometric -0.003749 0.0213 0.222 -0.1769 0.04717 


Velocity comparisons with the Navier-Stokes predictions were also 
encouraging. The velocity predicted by the boundary- layer equations 
compares well with the Navier-Stokes solution of Macagno and Hung [5] 
(Fig. 11) for Re = 60 and a 1:2 pipe expansion. Even in the region of 
reversed flow the agreement is good. Figure 12 shows the centerline 
velocities for axisymmetric flow for three different expansion ratios. 
The boundary- layer solution used for this plot is again that for just 
one sweep down the channel since global iteration did not affect this 
parameter. The agreement between the boundary- layer solution and the 
Navier-Stokes solution of Macagno and Hung [5] for the 1:2 expansion 
ratio is almost perfect. Pollard's [ 12 j predictions do not compare as 
well . 

The constant property boundary- layer solution, which is independent 
of Reynolds number (see Section II. B), can be thought of as the 
asymptotic limit to the solution of the Navier-Stokes equations as the 
Reynolds number becomes large. Figure 13 shows the Reynolds number at 
which the boundary- layer equations can be used in place of the Navier- 
Stokes equations for symmetric expansions with d/D - 1/2. For Reynolds 





BOUNDARY-LAYER, d/D=0.3 
POLLARD, Re=250, d/D=0.3 
POLLARD, Re=80, q/D=0.3 
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Figure 12. Comparison of centerline velocity for axisymmetric expansions predicted by the 
boundary- layer equations and the Navier-Stokes equations [5,12] 
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numbers below 100, the boundary- layer predictions of T^. are greatly in 

error. This trend is consistent with the T . comparison shown in Table 

mm 

5. For the case shown in Table 5 (Re = 37.3), the prediction of T *y 

the boundary -layer equations was 30% higher than that predicted by the 

Navier-Stokes equatitns; the X-location of Y . was 22% lower. Figure 8 
n min 

shows that the distance to the vortex center predicted by the boundary- 

layer equations is less than that predicted by the Navier-Stokes 

equations and measured experimentally by what appears to be a constant 

amount for all Reynolds numbers. The relative error becomes less at 

larger Reynolds numbers. Global iterative sweeps down the channel 

affected Y . by less than l v . 
min 

Figure 14 compares c^ predicted by the boundary-! ayer method and 
that predicted by Chiu [142]. Chiu used a partially-parabolized Navier- 
Stokes (PPNS) model that neglected the streamwise diffusion terms but 
included the elliptic effects of pressure and convection. The geometry 
was that of a 1:3 two-dimensional expansion; the Reynolds number was 39. 

Figure 14 shows that global iteration has little effect on the 
reattachment length (where c^ = 0.0). The first sweep down the channel 
using the boundary- layer equations predicts smaller absolute minimum 
values of c f in the reversed flow region than the PPNS equations. 

Global iteration produces the same minimum value of c^ as the PPNS 
predictions, but this value occurs at a smaller X. Figure 14 hints that 
the boundav- layer predictions tend to "squeeze" the region of flow 
reversal closer to the step. Hence, the c^ curve reaches a minimum 
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Figure 14. Comparison of friction coefficient predictions between the boundary-layer equations and 
the PPNS equations [142] for a 1:3 planar expansion 



132 


sooner than the PPNS model, and the distance to the vortex center is 
underpredicted as is shown by Fig. 8. 

Figure 14 also compares central differencing and upwind 
differencing of the X-derivative in the reversed flow region. A few 
test calculations with central differencing of the X-derivative in the 
X-convective term were compared with the local upwind differencing that 
was usually used when FLARE was not in effect. Figure 14 shows that the 
predictions are very similar for the two types of differencing. 

Figure 15 shows a comparison of c^ predicted by the boundary- layer 
equations and c^ predicted by Pollard [45] using the Navier-Stokes 
equations for a 1:2 pipe expansion. For less severe pipe expansions, as 
for the 1:2 expansion case, global iteration has very little effect, as 
is shown by Fig. 15. It should be noted that Pollard's prediction for 
the Re = 250 case shown in Fig. 15 overshoots the known fully-developed 
c^ Reynolds number product of 2.0 by 11%. Since Pollard provides very 
few details of his computational procedure other than that it is similar 
to the SIMPLE method of Patankar and Spalding [50] , it is not clear why 
c^ exceeds the fully-developed values in some cases. 

Figure 16 shows a plot of c ^ similar to Fig. 15 for d/D =0.7 and 
Re = 250 which is a relatively mild expansion and a high Reynolds 
number. For a Reynolds number this large, the agreement was expected to 
have been better. Pollard's c^ predictions do not support the 
supposition that the boundary- layer solution tends to '’push" the 
recirculation region closer to the step as does Chiu's [142] and Macagno 
and Hung ' s [ 5 ] . 
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Figure 15. Friction coefficient predicted by the boundary- layer equations with and without global 
iteration compared with the Navier-Stokes solution [45] for a 1:2 pipe expansion 
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Figure 16. Friction coefficient predicted by the boundary- layer equations compared with the 
Navier-Stokes solution (45] for a pipe expansion with d/I) =0.7 
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Figures 17 and 18 show predicted values of the skin-friction 
coefficient for different values of h/d for planar and axisymmetric 
expansions respectively. Both figures show results obtained using the 
"once-through" method with the FLARE approximation. When using these 
curves, one should realize that they many not be accurate in the 
recirculation region for low Reynolds numbers (less than 100). 

Figure 19 shows the effect of global iteration on the pressure 
gradient for a 1:3 two-dimensional expansion. The pressure gradient 
appears to be the parameter most affected by global iteration. The 
predicted dP/dX was very similar whether local upwind or central 
differencing of the X-derivative in the X-convective term was used when 
flow reversal was present. 

One interesting prediction of the boundary- layer equations using 
FLARE was the existence of a small secondary eddy in the corner formed 
by the wall and the step. For d/D = 0.5, this eddy was less than 
l/220th the length of the primary eddy. This eddy was discovered while 
using an extremely small AX for mesh refinement studies. The second 
eddy rotates in a direction opposite that of the large one. The flow 
situation very near the corner is similar to Stokes flow across the top 
of a wedge cut in a wall. The solution to this Stokes flow is a "stack" 
of eddies in the wedge, decreasing in size and intensity as one moves 
down in the wedge [143]. For runs employing global iteration, a grid 
fine enough to predict the secondary eddy was not used due to expense. 
(AX would had to have been less that l/400th the distance to 
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Figure 18. Friction coefficient predicted by the boundary- layer equations using FLARE for 
different expansion ratios for axisymmetric expansions 
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igure 19. Effect of global iteration on the dimensionless pressure gradient, d V (Uu . )dp/dx , for 
1:3 planar expansion (X=y./(d Re)) 1 
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reattachraent to place one marching station beyond the step In the 
secondary eddy for d/D = 0.5.) Those solving the Navier -Stoker 
equations usually do not use a grid fine enough to resolve the secondary 
corner eddy. 

c. Initial conditions Modifying the initial conditions along 
the face of the step as proposed by Acrivos and Schrader [27] to take 
into account the "collision velocity" at the face of the step as was 
discussed in Section II.A.4.b was found to be of minor importance. The 
algorithm to determine the necessary velocity at the step face [144] 
predicted nonzero velocities for moderately fine grids but predicted 
zero velocities as the grid was refined. To test i.he effect a 
"collision velocity" could have had on the solution, the effect of the 
following sinusoidal velocity along the face of the step was studied: 

U(0,y) = -A sin(~p), 0<y<h (4.2) 

where A is an arbitrary amplitude. The predictions of the boundary- 

layer equations with global iteration for flow through a 1:2 planar 

expansion was used as a test case. An extreme case A = C.15, which 

caused velocities at the step to be greater than those normally in the 

recirculation region, predicted the reattachment length, location of 

T . , and value of t , to within 2% of the A = 0 case. (For a 1:2 
man min 

expansion, the maximum speed of the fluid in the recirculation region 
was 0.12 in the negative X-direction. ) Figuie 20 shows that the 
difference between the c^ predicted by the A = 0 case and the A = 0.15 



Figure 20 
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case is not significant. Since a "collision velocity" was not predicted 
and would have had very little effect if it were predicted, a zero face 
velocity was used for all other calculations of the present study. 

d. Primitive variable results The primitive variable 
formulation was used to test methods for determining the pressure 
gradient in the channel. The most efficient method is similar to that 
used for the u-<j> variable formulation where the 2x2 me "rices below the 
diagonal are eliminated. Then, Eq. (3.43) with V^j = 0.0 is used to 
algebraically give an expression for & in terms of known coefficients. 
The results using this method were the same as those predicted by the 
U-f variable formulation. For internal flows, the U-V form is in fact 
somewhat easier to program. 

If one desires to use one of the available general block solvers, 
the channel mass flow constraint cannot be included in the algorithm 
before the back-substitution step without adding an extra equation used 
to determine the pressure gradient to each block of equations as was 
done by Cebeci [ 145 ] . If the extra equation is not added, an iterative 
procedure as discussed in Section III.B.2.b of the p^^vious chapter in 
required. The three different variations of the secant method 
investigated in this study had different levels of success. Variation 1 
was by far the superior algorithm. A planar expansion with a 1:2 
expansion provided a test case. 

2 ) Variation 2 As described in Section III.C.2.b, 
variation 1 nested the pressure secant loop in the Newton linearization 
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loop. This ensured that an accurate value of the dimensionless pressure 

gradient, (1, was obtained before making the next Newton linearization 

iteration. The secant iteration continued until the nerly predicted & 

had a relative change from the old value less than sc .ie pnsc-ibed 

limit. The relative difference between the first and second secant 

-4 '6 

predictions for & was less than 5.0x10 and usually lcsi than 5.0x10 
(The first secant iteration is the first prediction after the two 
calculations were done using the guessed values of &.) A relative 
change of 5.0x10 ^ was used as the criterion for convergence. Since 
only one secant itiration was usually sufficient, it appears that there 
is a nearly linear relationship between the pressure gradient error and 
the error in V at the centerline. 

A comparison of the results using variation 1 and the U-Y variable 
solution showed the velocities, pressure, and c^ agreeing to four 
significant digits. In general, the results were within 0.5% of the U-f 
scheme. There were no oscillations of c^ as predicted by Kwon and 
Pletcher [39] when the continuity and momentum equations were not 
coupled. The execution time war three times slower than the U-f scheme. 
This ratio seems reasonable since the pressure secant method had to 
solve the system of equations three times in order to find the correct 
&. 

The lumber of Newton linearization iterations was counted for each 
step down the channel. Generally, only one update of the provisional 
coefficients was necessary. For the smallest AX used, the relative 
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change between iterations was less than 0.0005 on the second calculation 
of the velocities at a given X-position. 

2) Variation 2 The Newton linearization loop inside the 
pressure loop caused the solution to converge to an erroneous solution 
and then diverge when predicting the expansion flow. It makes little 
sense to use Newton linearization iteration when the initial pressure 
derivative is inaccurate. The failure of variation 2 shows the 
importance of leaving the pressure derivative as a variable and solving 
for it along with the velocities at each X-station. Variation 2 
predicted the correct velocity profile when modeling the inlet flow 
between two parallel plates which involved no separation. 

3) Variation 3 Removing the Newton linearization loop 
from variation 1 t large oscillations for the first few X-stations 
beyond the step. Figure 21 compares the predictions of c^ for variation 
1 and variation 3 for a Reynolds number of 50. 

2. Heat transfer 

The results of this section are tor variable property flow unless 
specified otherwise. Both specified heat flux and specified wall 
temperature boundary condition predictions will be discussed. Global 
iteration was important when predicting the temperature field. 

a. Comparison with o ther data sets As mentioned in the 
Literature Review, there are no experimental laminar heat transfer data 
sets for symmetric rapid expansions. For this reason, comparisons were 
made with pipe entry flows to check the program. 
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Figure 22 compares the Nusselt number predictions of the present 
study with 5 constant property eigen- function solution [127J for air 
f owing in a pipe. The temperature profile is uniform at the inlet and 
the velocity profile is fully-developed. The boundary condition at the 
wall is given by a cor st ant wall temperature. 

Figure 23 shows the predicted bulk temperature, T^, for variable 
property inlet flow in a pipe with a constant specified heat flux at the 
wa~l. The velocity and temperature profiles at the inlet are both 
uniform. Figure 23 compares the finite-difference solution of the 
boundary -layer equations by Bankston and McEligct [146] with the present 
predictions using the ASHRAE property erpressions given in Appendix A 
and the power-law property expressions used by Bankston and McEligot. 

The the wall heat flux was 


V d/2) 

k.T. 


10 


1 x 

where k^ and T^ are respectively the thermal conductivity at the inlet 
and the absolute temperature at the inlet. 

b. Predictions Flow through a 1:2 pipe rapid expansion was 
predicted for a constant wall temperature and a constant wall heat flux. 
The inlet diameter in both cases was 50mm; the inlet temperature profile 
was uniferm at 10°C. 


11 Constant heat flux 


For the specified heat flux 


boundary condition predictions. Re = 200 based on inlet conditions, and 
= 25 V/m 2 . The constant heat flux prediction converged much slower 


F * ' ! t e 0 i f f ereno* 
E i gen F unot i on 



Figure 22. Predicted Nusselt number for inlet pipe flow having an initial fully developed v 
profile and uniform temperature profile compared with an analytical solution for 
constant wall temperature boundary condition [127] (Pr =0.7) 


Bankston & McEligot 
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than the specified wall temperature prediction. For the specified heat 
flux, the temperature field boundary conditions are specified by the 
slope everywhere except at the inlet above the step. Since the actual 
value of the temperature is "tied" down only for a small portion of the 
boundary, the solution for the specified heat flux boundary condition 
converges very slowly. Figure 24 shows the relative change of the wall 
temperature at three different x-values with respect to global 
iteration. The momentum and continuity equations were not solved each 
time the solution of the energy equation was obtained due to the 
computational expense. A hydrodynamic solution was computed on energy 
global iteration numbers 1, 3, 10, 30, and 45-50. The large "burp" in 
Fig. 24 at iteration 30 is caused by the u-rf/ values being recomputed. 
Figure 25 shows the effect of global iteration on the Nusselt number. 
There is very little difference between the 30th and 50th iterations. 
This indicates that 30 iterations are sufficient for convergence. 

Figure 26 shows the wall temperature and the bulk temperature computed 
by Eq. (2.£1) and Eq. (2.82). Ideally, the two bulk temperature curves 
should exactly correspond. 

2) Constant wall temperature The predictions for the 
const v , it wall temperature boundary condition were obtained for Re = 1000 
based on inlet conditions and T^ = 100°C. Only six global iterations 
were required for convergence of the specified wall temperature 
predictions. By specifying the temperature along the wall as opposed to 
the slope of the temperature as for the specified wall heat flux, the 



149 





Figure 24. Relative change of the wall temperature at three different channel locations with 
global iteration for a 1:2 pipe expansion 
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Figure 26. Wall temperature and bulk temperature for a 1:2 pipe expansion v th a heat flux 

boundary condition (Re^=200, q^=25 W/m 2 , Pr=0.7, inlet T^=10°C, x is dimensional) 
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solution responds to the boundary conditions much quicker. Figure 27 
shows the Nusselt number for the first and sixuh global iterations for a 
constant wall temperature boundary condition. Figure 28 shows the 
Nusselt number for three different mesh sizes after convergence. 

B. Turbulent 

This section describes the turbulent flow results. The differences 
between the predictions using the low-Reynolds -number k-e aquations and 
the high -Rev:. Ids -number k-E equations with the different near-wall 
models will be presented. 

The turbulent calculations were done using a Perkin Elmer 3240 
computer. An average run using a 121x120 grid required 3.5 minutes of 
CPU time for the first iteration and 1.2 minutes for subsequent 
iterations. The reason for the different CPU time requirements is due 
to the linearization procedure used after the first global iteration for 
the turbulent flow calculations. For the turbulent flow calculations, 
the Newton linearization was dropped after the first global iteration 
and the nonlinear terms were linearized using values from the previous 
global iteration. 

1. Hydrodynamic constant temperature 

a. Fully develope d pipe flow The hydrodynamic solution 
procedure was first tested by comparing with the well -documented fully 
developed pipe flow. Air at 17°C flowing in a pipe of internal radius 
of 0.1235m with a Reynolds number of 41,680 was used as a test case. 







Figure 23 Effect of grid refinement on the Nusselt number for a 1:2 pipe expansion with 
specified wall temperature (Re. =1000, T =100°C, Pr=0.7, inlet T =10°C) 



Uniform inlet conditions for u, k. and e were specified; the inlet 
condition for 0 was obtained by integrating the u profile. There were 
80 unequal grid divisions in y so that at least two grid points above 
the wall were in the viscous sublayer. 

Figure 29 shows the computer predictions of the velocity compared 
with Eqs. (2.60) and (267) in law-of -the-wall coordinates. Three 
different methods of spicifying the turbulent length scale were used in 
conjunction with the low-Reynolds -number k-equation: (1) an algebraic 

mixing length given by Eq. (2.56), (2) the high-Reynolds -number 
r -equation with the a boundary condition applied at y + = 30 and Eqs. 
(2.56) and (2.75) for y + < 30, and (3) the low-Reynolds -number 
e-equation solved throughout the flow field. Figure 29 shows that for 
fully developed pipe flow, all three methods give similar velocity 
predictions . 

Figures 30 and 31 show that the three different methods of 
specifying the near-wall mixing length predict widely varying values of 
k and Reynolds stress near the wall. The turbulent kinetic energy 
predicted by the low-Reynolds -number £ -equation is much closer to 
experimental values near the wall (Fig. 30). However, the Reynolds 
stress predictions using the low-Reynolds -number e -equation do not 
correspond to the experimentally measured values as well as those using 
the high-Reynolds -number £ -equation. 

For separated flow, the first global swe ip with the low-Reynolds - 
number model predicted reattachment lengths that were more than twice as 



Fig 7 



igure JO. fully developed pipe flow turbulent kinetic energy profiles for thru different lem ; 
scale models used with a low-Reynolds-nunber k-equation 
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long as measurements. It also proved to be unstable when iterating 
globally. Hwever, due to the good correlation between experimental 
measurements of k and those predicted by the low-Reynolds -number 
turbulence model, this model was used to provide the fully developed 
inlet values for u, if, k, and e for the step flow predictions using the 
hign-Reynolds -number k-e models. 

b. Rapid expans ion flow The laser anemometer measurements of 
flow through a symmetric planar expansion obtained by Smyth [64] were 
used as the test case for the constant temperature hydrodynamic 
predictions since he gives the most extensive set of data to date for 
symmetric expansion flows. For Tmyth's measur- ments, water flowed 
through a 1:1.5 planar expansion at a Reynolds number of 20,140 based on 
the inlet plate spacing and the inlet average velocity, u^. The flow 
was fully developed at the step. The initial conditions for the 
computer predictions were obtained by solving the channel flow case 
upstream of the step using the low -Reynolds -number k-e equations. An 
38x81 computational grid was used for th. computer predictions compared 
with the measurements of Smyth. 

Of the three different near-wall models mentioned in Section 
III.B.6.b, the maximum-sheai -stress model gave the best predictions. 

The variable-A + near-wall model gave generally better predic. ions than 
the constant-A + model in which the boundary conditions for k and t were 
given at y + = 30. Only the results of the maximum-shear-stress and 
variable-A + models will be presented for the constant temperature 


turbulent flow calculations. 
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All of the predictions for the flow except very near the wall were 
similar regardless of the near-wall model that was used. Watkins and 
Gooray [115] noticed similar behavior for their predictions using the 
Navier-Stokes equations. When they modified the wall functior used with 
the high-Reynolds -number turbulence model, the reattachment length and 
near-wall solution was affected but the predictions away from the vail 
were not. Since the predictions away from the wall are similar for all 
the near-wall models, they will be discussed below with the results of 
the maximum-shear-stress near-wall model. 

The predicted reattachment lengths were longer than that measured 
experimentally. The experimental reattachment po’it, deduced from the 
velocity profiles provided by Smyth [54], was 

1.2 < £ /r < 1.5 
r o 

The predicted re^t-techmene lengths according to iteration number and 
near-wall turbulence model were 

£ /r =2.2 variable-A + , first iteration 
r o 

1.8 variable-/,*, fiftieth iteration 
2.2 maximum-Reynolds -stress, first iteration 
2.1 me .mum -Reynolds -stress, fiftieth iteration 
Again, the reattachment length is a strong function of the near-wall 
turbulence model. 

The two main short comings of the variable-A model were slow 
convergence with respect to global iteration and the prediction of small 
unphysical irregularities near the wall. For the 88x81 grid used for 
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the comparisons with Smytl , 300 iterations wej e required for convergence 
of the variable-A + model. Figure 32 shows dp/dx at an x position 
approximately half-way to reattachment with respect to global iteration. 
Figure 32 shows that the computer predictions were not highly stable for 
the first 250 iterations even though all the variables were underrelaxed 
by 0.25 (the allowable change between global iterations was only 1/4 the 
predicted change). Due to the high number of iterations, the 
computational cost of solving the boundary- layer equations with this 
near-wall model is approaching the cost of solving the Navier-Stokes 
equations using coarse grids and wall functions near the wall. 

Figures 33 through 35 show the unphysical irregularities near the 
ca 1 for the variable-A + case. These irregularities are very noticeable 
in che velocity profile for x/r Q = 0.8 in Fig. 33, and the velocity and 
Reynolds stress profiles for x/r Q = 1.2 and 2.0 in Figs. 34 and 35. 
Figures 33 through 35 are the predictions for the converged solution 
after 600 iterations. 

For the predictions of the present study, the same mass flow down 
the channel was used as was reported in Smyth [64]. However, Figs. 33 
ard 34 suggest that for some experimental profiles (especially the inlet 
profile), more fluid is flowing down the channel than the predictions 
indicate. The experimental profiles were integrated using a trapezoidal 
rule find ; he mass flow rate in the channel. The mass flow measured 
experimentally was found t o v a .ry from 1% below to 13% above that 
r '.ported by Smyth, If anything, the trapezoidti rule integration would 
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tend to unde vpvedict the mass flow rate. Since Srnyth gives the accuracy 
of his u/^ measurements as ± 0.05, this author is uncertain of the 
cause of the discrepancy. 

The wall functions based on the maximum Reynolds stress with the 
boundary conditions applied at a y position corresponding to a y + of 30 
to 50 for a fully developed turbulent profile gave the overall best 
results. This method converged much faster than the variable-A + method. 
Only 50 global iterations were required for convergence with ar 
underrelaxation factor of 0.5. The unphysical irregularities in the 
velocity profiles were not predicted with the maximum-shear-stress near- 
wall model. However, irregularities ir> the Reynolds stress profiles 
were still predicted near the wall. 

Figures 36 through 39 compare the velocity profits for the first 
and last (fiftieth) global iterations with the experimental measurements 
of Smyth [64]. The predictions for the last global iteration (Figs. 38 
and 39 ^ are generally better than the predictions for the first (Figs. 

36 and 37) but not remarkably so. Other than in the separated region 
itself, global iteration has no noticeable effect cn the velocity 
profiles. This indicates that the flow outside the separation bubble is 
not affected by the way the x-co ivective term is approximated. The 
re ison for the different predictions for the first and last global 
iterations in the separation bubble is due to convection in the negative 
x-directior. in this region. The velocity predictions on the first and 
the last global iterations in the recovery region downstream of 





Figure 37. Compariso; 

model and 







0 '0 



profile predicts 
d by Smyth [64] , 








171 


reattachment are veiy similar to each other and are not in good 
agreement with experimental measurements. 

Figures 40 and 41 compare the predictions of the turbulent kinetic 
energy on the first and last global iterations with the experimental 
measurements of Smyth. The predicted turbulent kinetic energy varied 
little between the first and last global iterations except near the 
step. The maximum k predicted numerically is less than that measured 
experimentally by Smyth. The k values also drop off too rapidly as the 
channel centerline is approached. This last point is not too important 
because near the centerline, 3u/3y is small so any error in the 
turbulfnt viscosity will not have a very large effect. 

Global iteration has some effect on the Reynolds stress as is shown 
by Figs. 42 and 43. On the last iteration, the Reynolds stresses were 
observed to be higher in the recirculation region (x/r Q = 0.4) than on 
the first global iteration. The distance from the wall to the point of 
the peak Reynolds stress does not decrease as reattachment is approached 
and then increase after reattachment has been passed as experiments have 
shown should be the case [1]. The maximum Reynolds stress is also too 
large after reattachment in the redeveloping boundary layer. 

2. Heat transfer 

The turbulent heat transfer results are discussed in this section. 
For the predictions, the variable property relations for air given in 
Appendix A were used. The computational field had a 120x121 grid for 
the heat transfer predictions. The predictions of the present study 
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Figure 42. Comparison of Reynolds stress profile predictions using the raaximum-shear-stress near- 
wall model and those measured by Smyth [64], first iteration (d/D=2/3, Re -20140) 
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will be compared with the Nusselt number measurements of Zemanick and 
Dougall [84] and Baughn et al. [90], and the numerical predictions of 
Johnson and Launder [113] and Watkins and Gooray [115]. 

To verify the computer program, the fully developed temperature 
profile in law-of-the-wall coordinates was calculated for pipe flow. 

The dimensionless temperature measured experimentally in a fully 
developed flat plate turbulent boundary layer [127] is compared with the 
predictions for fully developed turbulent pipe flow in Fig. 44. It 
would have been better to compare with experimental measurements of 
fully developed pipe flow but none were readily available. The velocity 
law-of-the-wall profiles for boundary- layer flow and fully developed 
pipe flow are very similar. For this reason, the temperature law-of- 
wall profiles for the two flows should also be very similar in the 
logarithmic region, although one might expect minor differences. The 
dimensionless temperature, T + is given by 

T + ■ <>y T „ - T >v% 

The agreement between the two curves is acceptable. 

The variable turbulent Prandtl number developed by Watkins and 
Gooray [115] reduces to a value near 0.2 for fully developed turbulent 
pipe flow (Fq. 2.81). Thus, one would expect the Pr^ expression of 
Watkins and Gooray to greatly increase the heat transfer predictions 
over those predicted using a constant Pr of 0.9. When the Pr 
expression of Watkins and Gooray was used in the computer program of the 
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present study, the heat transfer predictions in the separated flow were 
generally more than 100% higher than the experimental measurements. 

Since predictions using the variable Pr^ expression given by Eq. (2.81) 
predicted extremely high heat transfer, the results using a variable Pr^ 
will not be discussed further. The remaining heat transfer predictions 
are for Pr^ equal to 0.9. Several other recent numerical predictions 
successfully used a constant Pr of 0.9 when predicting the heat 
transfer in a pipe expansion using the Navier-Stokes equations [89, 113, 
119] . 

The maximum-shear-stress inner viscosity model giver, by Eqs . 

(2.65), (2.76), and (2.77) not only gave the best hydrodynamic 
predictions but also gave the best predictions for heat transfer. The 
Nusselt number predictions obtained using the maximum-shear-stress 
iimer-model and the constant-A + model are compared to experimental 
measurements [84, 90] in Fig. 45 for d/D of 0.8 and Re near 20,000. 
Figure 45 shows that global iteration improves the heat transfer 
predictions in the recirculating flow. The overall agreement is quite 
good considering the difficulty in predicting turbulent heat tiansfer in 
separated flew. 

Figure 46 is a comparison similar to Fig. 45 but for a more extreme 
expansion, d/D - 0.53, and a lower Reynolds number. Re = 10950. The 
predictions of Watkins and Gooray [115] and Johnson and Launder [113] 
for the same flew case are also included. For this large expansion, it 
was very difficult to specify initial values for k and t along the face 
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igure 45. Comparison of the predicted Nusselt number using the Prandtl-mixing-length (A -25) and 
maximum-shear-stress near-wall models with experimental measurements [84,90] tor an 
axisymmetric expansion (d/D = 0.8, Re = 20161, Nu based on D) 



Measurements * Baughn et al. , Re 0 = 10954 

+ Zemanick & Dougal 1, Ra o -9620 

Predictions — Present, 1st Iter. 

present, 50th Iter. 



Figure 46. Comparison of the predicted Nusselt number using the maximum-shear-stress near-wall 
model with Navier-Stokes predictions [113,115] and with experimental measurements 
[84,90] for an axisymmetric expansion (d/D = 0.533, Re = 10954, Nu based on D) 
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of the step so that the prediction algorithm was stable for the first 
sweep down the channel. For the algorithm to be stable, the variable 
expression for c^ (Eq. 2.80) had to be discarded in favor of a constant 
value of 0.09. Figure 46 shows that the heat transfer is overpredicted 
in the recirculating region and in the initial redeveloping boundary 
layer but is then underpredictad as the boundary layer develops. Again, 
as was indicated by Fig. 45, global iteration greatly affects the heat 
transfer predictions in the reversed flow region. 

The predicted x-pos l. ion of Nu , x/h = 5.4 for the fiftieth 
1 r max 

global iteration, wa. : near the value measured by Baughn et al . [90] for 
d/D - 0.8 (F’g 45). The predicted reattachment point was l / h = 6.0, 
slightly downs ti earn of the point of maximum Nu. For d/D = 0.553, the 
p.adit ri rt. attachment point was upstream of the predicted point of 
maximum Nu. The predicted reattachment length and point of maximum Nu 
were 9 2 and 10.1, respectively. This predicted point of maximum Nu was 
near the value measured by Baughn et al. but downstream of the other 


studies shown in Fig. 46. 
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V. CONCLUSIONS 
A. Laminar 

From the results of the laminar predictions mentioned in the 
previous section, one can draw the following conclusions. 

1. The distance to flow reattachment and velocities outside of 
the trapped eddy are very well -predicted by the boundary- 
layer equations for Reynolds numbers above 20 and expansion 
ratios below 1:3. 

2. The eddy structure is not wel 1 -predicted for low Reynolds 
numbers (Re < 200 for a 1:2 pipe and planar expansion). This 
is shown by the poor predictions of the magnitude and 
position of ^ m ^ n - For planar expansions, when the Reynolds 
number is approaching the point where the eddy structure can 
be predicted by the boundary- layer equations, experiments 
have shown the flow tends to be either asymmetric or unsteady 
[ 10 ]. 

3. The flow can be qualitatively divided into two regions. An 
elliptic region made up of the trapped eddy and a parabolic 
region including the rest of the flow field. The eddy 
asserts only a weak influence on the rest or the typically 
parabolic flow. 

4. Global iteration over the flow field using the boundary- layer 
equations does not significantly change the hydrodynamic 
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constant temperature predictions and offers little 
improvement over the "once through" method using FLARE. 
However, global iteration is important when predicting the 
heat transfer for an insulated step face, as was the case for 
the present study. When the heat flux is specified at the 
wall, the convergence of the global iteration to determine 
the temperature field is approximately five .imes slower than 
when the wall temperature is specified as a boundary 
condition . 

5. A zero velocity initial condition on the face of the step is 
adequate. A nonzero velocity on the step face was not 
predicted by the algorithm of Acrivos and Schrader [27] as 
the grid was refined. A small nonzero velocity on the face 
of the step did not significantly affect the solution. 

6. The primitive variable formulation (U-V) predictions were the 
same as the U-t formulation predictions. Thus, the "wiggles" 
previously noted with primitive variable solutions [39], were 
due to uncoupling of the continuity and momentum equations 
and not the choice of variables. 

B. Turbulent 

This section describes the conclusions that can be drawn from the 
turbulent hydrodynamic and heat transfer predictions. 
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1. Although the low-Reynolds -number k-e turbulence model of 
Chien [135) gave excellent results for attached flow in 
channels and pipes, it greatly overpredicted the reattachment 
length wh*n used for rapid expansion geometries . The low- 
Reynolds -number model was also unstable when iterating 
globally. The high-Reynolds-number k-c turbulence model with 
near-wall turbulence models gave better predictions. 

2. The near-wall turbulence model exerted a very weak influence 
on the flew field away from the wall. It had a very strong 
influence on the near-wall flow and the point of reattachment 
as was reported by Watkins and Gooray [115]. The near-wall 
model had a primary role in determining the heat transfer 
rate. 

3. The near-wall maximum-shear-stress turbulent viscosity model 
based on the inner viscosity model of Johnson and King [121] 
gave better predictions than the near-wall models based on 
Prandtl's mixing length. The maximum-shear-stress model 
required less underrelaxation when iterating globally, 
required fewer global iterations to converge, and predicted 
the heat transfer much better than the Prandtl mixing length 
near-wall models. The maximum-shear-stress model did not 
predict small irregularities in the velocity profiles near 
the wall as did the variable-A + model. 
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4. A modification of the turbulent Prandtl number expression, 
Pr t , derived by Watkins and Gooray [115] when used in the 
computer code of the present study, predicted values of Pr 
between 0.2 and 0.3. This caused the heat transfer 
predictions for attached pipe flow to be greatly 
overpredicted. The near-wall models based on Prandtl' s 
mixing length underpredicted the heat transfer in the 
recirculating region when using Pr^ = 0.9. Using the 
expression for Pr^ of Watkins and Gooray with the near-wall 
models based on Prandtl' s mixing length for the separated 
flow may give reasonable heat transfer predictions in 
separated flow but does not appear to be applicable for flows 
with no separation. 

5. Global iteration, as for laminar flow, was not very 
important for the hydrodynamic predictions. However, it was 
even more important than for laminar flow when predicting the 
heat transfer. 

6. The computer algorithm was sensitive to the initial values of 
k and t specified at the step necessary to start the first 
sweep down the channel. The larger the expansion, the more 
critical the initial values of k and e became. 

7. As for laminar flow, the parabolic region oi the flow outside 
the separation bubble was not affected by global iteration. 
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8. The predictions of the peak turbulent kinetic energy were 
less than those measured by Smyth [64] . 

9. The overprediction of the Reynolds stress in the 
redevelopment region caused the underprediction of u in this 
region. The predicted peak value of -u'v' did not dip toward 
the wall as reattachment was approached and then move away 
from the wall as it was passed as indicated by experiments 
[ 1 ] - 

10. For the less extreme of the two pipe expansions examined, the 
heat transfer was well-predicted in the separated region and 
slightly overpredicted as the flow developed. For the more 
extreme expansion, the heat transfer was overpredicted in the 
separated region and near reattachment. 

11. The boundary- layer equation method of this study provided a 
relatively inexpensive way to evaluate turbulence models 
applicable to separated flow that occurs with such devices as 
turbines, heat exchangers, airfoils, and ramjets. 

C. Recommendations for Future Research 

The present study indicated that the near-wall model derived from 
the inner viscosity model of Johnson and King [121] improved the heat 
transfer predictions. It may be desirable to study the use of this 
near-wall model based on the maximum Reynolds stress with the Navier- 
Stokes equations to predict the heat transfer in separated flow. Since 
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the correct prediction of the heat transfer rate is so sensitive to the 
near-wall model, heat transfer predictions should not be ignored when 
developing turbulent models valid in reversed flow. 

The use of wall functions for velocity for the near-wall flow, in 
contrast to solving the momentum, continuity, and energy equations up to 
the wall as was done in the present study, is becoming widely accepted 
as the most economical and robust method of predicting turbulent flow 
with separation. However, these three PDE's are valid near the wall and 
might be used to provide insight to turbulence modeling of the wall 
dominated region. The fact that the wall is not approached with the 
wall-function methods points to the lack of the ability to model the 
turbulence in this region. More research should be conducted in solving 
the actual momentum, continuity, and energy partial-differential 
equations with the turbulence model input as either a turbulent 
viscosity or the Reynolds stresses themselves in the near-wall region. 
This would clearly show the shortcomings of the near-wall turbulence 
models and hopefully lead to more universal ones. It is understood that 
solving the momentum, continuity, and energy equations in the near-wall 
region may not be the most cost effective way to predict turbulent flew 
at the present time, but information from doing so may help in the 
development of more universal wall models. 
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VIII. APPENDIX A: VARIABLE PROPERTIES FOR AIR AND WATER 


The equations used to approximate the properties of air and water 
as a function of temperature were taken from Thermophysical Properties 
of Refrigerants { 145 ] . 


A. Air 


The density of air is given by the perfect gas law 

pCkg/m 3 ) = (A. 1) 

wher?. T is the absolute temperature, R is 287.0 J/(kg*K), and p is the 
pressure in Pascals. 

'The viscosity is approximated by 

U(10-‘Ns/m 2 ) = /f/(0. 671692 + 85.22974/T - 2111.475/T 2 

+ 106417/T 3 ) (A. 2) 


for 60 < T < 1000 K. Equation (A. 2) has a maximum deviation of ± 0.7% 
and an average deviation of ± 0.2%. 

The thermal conductivity is given by 


k(W/(m*K)) = A/ (A + B/T + C/T 2 + D/T 3 ) 


(A. 3) 


where A, B, C, and D have the following values 

Range (K) A B CD 


80-300 

385.859 

9.11440x10* 

-2.68667x10* 

5 . 52604x 10 7 

300-600 

328.052 

1.67320x10* 

-3 02953xl0 7 

3.05682x10’ 

600-1000 

539.544 

-3.32903x10* 

3.59756x10* 

-9.67202xl0 1# 
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The maximum error of Equation (A. 3) is ± 0.23%. 

The specific heat at one atmosphere is given by 


c p (kJ/(kg*K)) = A + BT + CT* + DT 3 


(A. 4) 


The constants in Eq. (A. 4) have the following values 


Range (K) A B CD 

90-260 1.03200 -1.22500x10"* 0.0 0.0 

260-610 1.04466 -3.15967x10"* 7.07909xl0" 7 -2.70340xl0" x# 

610-900 1.00205 -1.62983x10"* 5.69525xl0“ 7 -2.68081xlO" 10 

Equation (A. 4) is accurate to within ± 0.018% for T > 260 K and accurate 
to within ±0.8% for T < 260 K. 


B. Water 


Water is considered an incompressible fluid with p equal to 
995.6 kg/m 3 . 

The viscosity of water is highly temperature dependent. The 
recommended equation is 


y(10~ 3 Ns/m 2 ) 


(A + B/T + C/T 2 ) 
e 


(A. 5) 


the constants are 


Range (K) 

A 

B 

C 

273-350 

0.030185 

-2191.60 

6.38605x10* 

350-500 

-3.22950 

13.18574 

2.65531x10* 

500-620 

-8.77361 

5875. P" 

-1.28275x10* 


The maximum deviation of Eq. (A. 5) is ± 1.5% from measured values. 



204 


Tht thermal conductivity of water can be approximated by 


k(W/(m*K)> = A + BT + CT 2 + DT 3 


(A. 6) 


to within ± 0.19% with the following constants 

Range (K) A B C D 

273-400 -0.61694 7.17851xl0" 3 -1.16700x10“* 4.70358x10"* 

400-600 -0.14532 4.02217xl0" 3 -4.64993x10"* -4.89257xl0" l# 


The following equation gives tue specific heat of water for 237 < T 
< 450 K to within ± 0.13%. 


c (kJ/(kg*K)) = 17.6611 - 0.147914T + 6 . 086 J 9x lO"*! 2 
P 

-1. 11867xlO"*T 3 + 7 . 8G297x 10" 1 °T fc (A. 7) 
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IX. APPENDIX B: RESULTING COEFFICIENTS FROM THE DISCRETIZATION OF THE 

MOMENTUM AND CONTINUITY EQUATIONS (u-* VARIABLES) 

The coefficients in the continuity equation are 


e j = b j = l (p j 1 -4 r j-i 4y - ) 


The coefficients of Eq. (3.14) when the FLARE approximation is in use or 
when u is nonnegative (j/NJ) are 

B j ■ 

c j ■ + sfj^j-1 - • B V] 

j - zrh 1(2fl j - “J> • ■ ♦/>]♦ 


D. = 


— (^Tx M 


B 2 


8j '‘Ay* i j+i“j+i (Ay + +Ay_) M j r j + Ay? r j-$ M J-i^ 

E = — * — C^-*u - - Bu ) 

E j 6 j Ax_ ( Ay_ u j-l Ay + j+1 ;T 

H. = 1.0 
J 

S C Ay^ Ay + j 
0j = r (Ay + + Ay) 


(B.l) 
(B . 2) 


After the first global iteration, if u is negative the coefficients 
of Eq . (3 . 14) are 

a j ■ ri 7^ ( V*r 2) ■ zr>-ViYi + Yj«v*->|) 
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b j - dri&yrv - 4:K r j-i M j-i ■ Vi (4 v‘ y - 
c j = srJ’i'V 1 + • Sfr+i ' 6 V) 

d j 1 ^[ p j fu j i+2 • *V + Nj - *j lt2) ] + 

• (IyJ«y.) M j r j + ^! r j-i M j-i ) 

e = 

j 6jAx + Ay_ j-1 Ay + j+1 


- * > 


H. = 1.0 
J 


3 and 8 . are the same as in Eqs. (B.l) end (B.2). 
J 

For j=NJ, the coefficients of Eq. (3.14) are 


Arj = E NJ = °-° 
2 m+l 

®NJ = " Ayl M NJ-i 

c nj = aT p nj^nj )2 


,m+l 


n 1 i* , 2 u 

D., T * — p • ( 2u v . T - u., T ) + 7-5 M klT . 

NJ Ax j NJ NJ Ay; NJ-j 


NJ 


1.0 
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X. APPENDIX C: MODIFIED THOMAS ALGORITHM 
Upon reducing Eq. (3.20) to an upper triangular matrix, the 

£ it 

diagonal submatrices [D] . and right hand vector of knowns {0} are 

J J 



* •’-> 
a 

i 

F." 


* * 
H.X + C, 

* 

[D] = 

J 

; (O* = 

J J 

j 

* 

-1 . 

j 

G.X + d. 
J J 


Ej remains unchanged from the value given in Appendix B. The modified 
coefficients in terms of the coefficients of Appendix B are (j=2,NJ) 


= 2CB .A . , + D. 

J J-l J 

= »A._i(b. + .*.j) + e. 


(C . 2 ) 
(C.3) 

it 



= YB . (H , , + E. ,G. .) + H. 
J J-l J-l J-l J 


(C.4) 

* 



= YB.(C. . E. .d. ,) + C. 
J J-l J-l J-l J 


(C . 5) 

* * 

it 


= ♦ .J.P + ®j- 1 (bjEj_ 1 - 

Vi» 

(C.6) 

* it 

it 


- UCj.jCb. ♦ Vl ) + dj.^b.Ej.j - 

u j-i» 

(C . 7) 

* * 





(C.8) 

are unchanged by the transformation 

to upper triangular 

form 


so are not starred. 

The boundary conditions at the wall are used to determine the 
coefficients for j=l. The following is true for j=l, 




* r T i+l, r . i+1, 

D i E i “i r*, ° u 2 

. + 

* i+l i 

e x -lJUj LO Ojl^ 


* >v 

H l x+C l 


i+l i+l * * 

Since - 0, D^, E^, and are arbitrary, they are set 

equal to 1. A^, H^, C^, c^, G^, and d^ = 0 cause Eq. (C.9) to be 
satisfied. Now that the coefficients for j=l are specified, Eqs . (C.2) 
through (C.8) can be evaluated starting with j=2 and continuing until 
j=NJ. 

Solving the reduced continuity and momentum equations for j=NJ (E4. 
3.24) gives the following expressions for and X. 

i + 1 *■ * * * * * 

U NJ = f(E NJ G NJ +H NJ ) ^Nj“ G NJ C NJ +d NJ H NJ ]/(e NJ H Nj‘ G NJ D NJ ) 

* i+l * 

* = (e N.T u NJ ' *NJ ‘ d NJ )/G NJ 

Back substitution is now used to solve for the unknowns u^ + ^ and 

J J 

for 2<j<NJ-l with the following expressions 


[ (H*+E ,G . )x - 
J J J 


i+l * * * 

A.u. . , + C. + E.d.l (D.+E.e.) 
J J+l J J J J J J 


* i+l 
e .u . 

J J 
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XI. APPENDIX D: RESULTING COEFFICIENTS FROM THE DISCRETIZATION OF THE 

ENERGY, k, AND t EQUATIONS 

The coefficients of Eqs. ("*.26) and (3.31) using the FLARE 
approximation ard when u is positive are 





1 

y»+ 

L 

-f A 

0 .Av A: 

J ' 


1 

i 

* — C P - u 
Ax_ r j , 

_ 1 j 

i 

Ax_ 

cp J 

L j 

2 1 

r 

■b + 

it 

Ay_ : 

r 2 

Ay. 

Ay+ 


>,c 


r J+iVj+t + CA V^ r 0 ,j r ij 


0,d 

\ 


P = — 
*3 


2 | ~Ay p 1 p 2 r r a. r j 

j+iVj+i Ay + +Ay_ P 1 0 , j j Ayf. j~i 


where & and 0^ are the same as in Eqs. (B.l) and (B.2) and <f> takes the 
value of H, k or e depending oi the transport equation being solved. 

The source terms for the cases when 0 is k or e are given in 
Section III. B. 6.1. For the energy equation, 

s h,c * r(^ r j+4Vi[ c Vi , ’’ <u j )1 ] + 


S H,d = °-° 
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where L = as given in Table 4, 

After the first iteration when the FLARE approximation was not 
used, the following coefficients apply when u < 0: 


* h] 

b i = ^ + <s i+i > - 

1 i i+1 ,i+2 . _ 
c. = - t — p u. i . + S . 

J J J J 


. _ 1 r i i+1 B i+2 .i+1.1 t 

a . - - t — p.u. + 'sriit. ~4> . ) + 

J Ax+L J J e j J J J 


P 3 ' %.d 


P_, P_, S. and S. , are defined above. 

2. 3 y jC ? 

For j=NJ, the following set of coefficients apply 


“NJ = 0 0 
. _ 1 „m+l. 


NJ 


Ay: 


t2“ 




- 1 i i+l.i , o 

c. . t 7 — p.u. 0. + 

NJ Ax/j j r j #,c 

, _ 1 i i+1 . 1 ^nr+lp 

NJ Ax_ p j U j Ay^ 2 *,j-i S 0,d 

Equation (3.30) implements the wall boundary condition for the 
energy equation when a heat flux is specified. The coefficients are 
obtained from a three-point finite-difference approximation of 3T/3y and 
are given as follows: 

2Av,+Ay 

j _ 1 ~ _ . 

1 ^y 1 (Ay 1 +Ay 2 ) 


Ay t +Ay 2 
1 i yi Ay 2 


a. = 
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where 


Ay i 

Ay 2^ 4y i +ly 2* 



, 1. , i+l*2 

+ 2 a l u 2 > 


XCufV] 
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XII. APPENDIX E: RESULTING COEFFICIENTS FROM THE DISCRETIZATION OF THE 
MOMENTUM AND CONTINUITY EQUATIONS (U-V VARIABLES) 

The coefficients of Eq. (3.42) are 


A. = 


(V, - 7~) 


j AY + +AY_ AY + 


B ss - (V - -?— ) 

j AY + +AY_ k j AY/ 


= ^-C2U ; -U") + 


2 (-i- + -i-) 


j AX' j j' AY +AY/AY AY_ 

(*' A - 1 A A A A 

c j ■ sr ( V - s^r ( -Vi 9 j + YiY 


E. = 


j AY + +AY_ (U j+l ' °j-l ) 


H. = 1.0 
J 


. - 

b j " e j 2AX_ 


dj - b.cuj ♦ U^) 


G. = 0.0 
J 


(E-1) 
(E.2) 
CE.3) 
(E.4) 
(E .5) 
(E.6) 
(E. 7) 
(E.8) 
(E .9) 


Equations (E.l) through (E.9) are valid for j=l to NJ except that is 

twice the value given by Eq. (E.2). 

After the system of equations (Eq. 3.41) is reduced to upper 

* * 

triangular form, the elements of [D]^ and {CJ^ are 


* 

D = D. 
j J 


£a. .B. 
* J-l J 


•j = *i ' + e j-i> 
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V V 'p at - o^.p 


c > C J - *v c j-l • 


G% G. - - oJ.jIj.j, ♦ 

d ! ■ d j - i'VVi - -w* w > + 'vtf-i - V 1 Vi )l 


* * 

I = D. - E. ,e. , 

j-1 j-1 j-1 


The above equations give the modified coefficients at the j level 
as a function of those at the j-1 level. A. and E. are not affected by 

J j 

the change to upper triangular form. Those for j=l are found from the 

i+1 i+1 ****** 

boundary conditions =V^ =0. Taking A^, Hj , Cp G n , d^, Ej, 

* * 

and e^ all equal fc ziro and Dj equal to one satisfies the boundary 
conditions . 

After reducing the coefficient matrix to block upper triangular 
form, Eq. (3.43) with j=NJ is used to solve for U*** and $. The 
result is 


U NJ (d NJ H N J^NjSj* ' / (e NJ H NJ _D NJ G NJ ) 


6 (e NJ U NJ ’ d NJ )/G NJ 


After knowing P and back substitution can be used to find 


the rest of the Us and u s for j=NJ-l to j=2. 


,.i+l , ,,i+l 

U. , and V. are given 
J J 


by 


214 


u j +1 = a I_A j u j+i + (H j ‘ E j G j )e + c j * E j d *J 


V i+1 = lfe*A U i+1 4 . for * * * 

j 2‘ j jVl (0G j ‘ ®j (H j • E j G j))& + 


d ? ' V c * - E j<yi 


_ * it 

2 = D. - E. e . 
J J J 
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XIII. APPENDIX F: Y-GRID STRETCHING TRANSFORMATIONS 


The inverses of the general stretching transformations of Roberts 
[138] as cited in Anderson et al. [4] were used to transform a uniform 
grid spacing ((-space, A( = constant) to a nonuniform grid spacing (y- 
space). The x = 0 line was divided into two regions, the step face and 
the inlet region. In both of these regions, 0 < ( < 1. 

Along the face of the step, 0 5 y S h, the transformation was given 
by 



At the inlet, h < y S y Nj’ t * ie trans f orma t* on was 


y = 


r 

0+1 

- (0-1) 

FI) 

(l-c)' 

1 1 

mr 

+ 1 



+ h 


The values of NJ, o, and the number of grid points below the lip of 
the step (y = h) were adjusted until there was a smooth variation of Ay 
near y = h and enough grid points near the wall to resolve the laminar 
sublayer regions. The values of o typically ranged from 1.005 to 1.05. 


